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Preface 



The last thing one settles in writing a book is what one should put in first. 

Pascal's Penates, 19 

In 1902 Jacques Hadamard introduced the term well-posed problem. His 
definition, an abstraction from the known properties of the classical problems of 
mathematical physics, had three elements: 

Existence: the problem has a solution 

Uniqueness: the problem has only one solution 

Continuity: the solution is a continuous function of the data. 

Much of the research into theoretical physics and engineering before and 
after 1902 has concentrated on formulating problems, with properly chosen initial 
and/or boundary conditions, so that their solutions do have these characteristics: 
the problems are well posed. 

Over the years it began to be recognized that there were important and 
apparently sensible questions that could be asked that did not fall into the 
category of well-posed problems. They were eventually called ill-posed problems. 
Many of these problems looked like a classical problem except that the roles of 
known and unknown quantities had been reversed: the data, the known, were 
related to the outcome, the solution of a classical problem; while the unknowns 
were related to the data for the classical problem: they were thus called inverse 
problems, in contrast to the direct classical problems. (Later reflection suggested 
that the choice of which to be called direct and which to be called inverse was 
partly a historical accident.) For completeness, one should add that not all 
such inverse problems are ill-posed, and not all ill-posed problems are inverse 
problems! This book is about inverse problems in vibration, and many of these 
problems are ill-posed because they fail to satisfy one or more of Hadamard's 
criteria: they may not have a solution at all, unless the data are properly chosen; 
they may have many solutions; the solution may not be a continuous function 
of the data, in particular, as the data are varied by small amounts, it can leave 
the feasible region in which there is one or more solutions, and enter the region 
where there is no solution. 
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Classical vibration theory is concerned, in large part, with the infinitesimal 
undamped free vibration of various discrete or continuous bodies. This book is 
concerned only with such classical \ibration theory. One of the basic problems 
in this theory is the determination of the natural frequencies (eigenfrequeneies 
or simply eigenvalues) and normal modes of the vibrating body. A body that is 
modelled as a discrete system of rigid masses, rigid rods, massless springs, or as 
a finite element model (FEM) will be governed by an ordinary matrix differential 
equation in time t with constant coefficients. It will have a finite number of 
eigenvalues, and the normal modes will appear as vectors, called eigenvectors. A 
body that is modelled as a continuous system will be governed by a set of partial 
differential equations in time and one or more spatial variables. It will have an 
infinity of eigenvalues, and the normal modes will be f mictions, eigenfunctions, 
of the space variables. 

In the context of classical theory, inverse problems are concerned with the 
construction of a model of a given type, i.e., a mass-spring system, a string, 
etc., that has given eigenvalues and/or eigenvectors or eigenfunctions, i.e., given 
spectral data. In general, if some such spectral data are given, there can be 
no system, a unique system, or many systems, having these properties. In the 
original, 1986, edition of this book, we were concerned exclusively with a stricter 
class of inverse problems, the a> called reconstruction problems. Here the data 
are such that there is only one vibrating system of the specified type which 
has the given spectral properties. In this new- edition we have widened the 
scope of our study to include inverse problems that do not fall under this strict 
classification. 

Before describing what the book is, we first say what it is not: it is not 
a book about computation. In Engineering, the almost universal approach to 
inverse problems is through least squares: find a system which minimizes the 
distance between the predicted and desired behaviours. While the early studies 
were examples of brute force, there is now an established and rigorous discipline 
governing such approaches, based on the work of Tikhonov, Morozov etc. See for 
example Kirsch (1996). We do not refer to any of this work in this book. Rather, 
we are concerned with basic analysis, qualitative properties, whether a problem 
has one or more solutions, etc. There are occasions when one method that we 
describe, that should theoretically lead to the construction of a solution, is found 
in practice to be ill-conditioned, and this has led to another, better behaved, 
procedure; in such a case we have presented both methods and discussed why 
one fails while the other succeeds; see for example Section 4.3. Because we 
are concerned with fundamental analysis, the range of physical systems that we 
can consider is relatively narrow; essentially it is confined to the basic elements 
of structures, rods, beams and membranes, and excludes structures composed 
of combinations of these elements. This restriction in scope is understandable; 
indeed, until the introduction of the finite element method and high-speed large- 
memory computing, the only direct vibration problems that could be solved were 
those involving those same structural elements in isolation. The study of inverse 
problems is at an earlier stage of evolution than that of direct problems. 
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The book fails into two parts: Chapters 1-9 are concerned with discrete 
systems, Chapter 10-14 with continuous systems. 

Matrix analysis is the language of discrete systems, and it is developed, as 
needed, in Chapters 1 and 3. Thus, Chapter 1 provides the bask definitions 
and introduces quadratic forms, minimax theorems, eigenvalues, etc. Chapter 2 
provides the basic physics of the vibrating systems that are analysed. Chapter 
3 lays out the classical analysis of Jacobi matrices, the matrices that appear in 
the simplest kinds of vibrating systems, in-line sequences of masses connected 
by springs. Chapter 4 concerns inverse problems for Jacobi matrkes. Chapter 
5 provides an introduction to more general discrete systems, and the language 
of graph theory that is needed to analyse them. 

Inverse problems in vibration are concerned with constructing a vibrating 
system of a particular type, e.g., a string, a beam, a membrane, that has speci- 
fied (behavioural) properties. The system so constructed must be realistic: its 
defining parameters, masses, lengths, stiffnesses, etc., must be positive ; Signs, 
poatfre and negative, lie at the heart of any deep discussion of inverse problems. 
Chapter 6, on Positivity, introduces the mathematics relating to different kinds 
of matrices: positive, totally positive, oscillatory, etc. This mathematks, due 
to Fekete, Perron. Gantmacher, Krein and others, was first applied to vibrating 
systems by Gantmacher and Krein in their classic Oscillation Matrices and Ker- 
nels and Small Oscillations of Mechanical Systems (1950), that has just recently 
(2002) been reprinted by the American Mathematical Society. 

Sometimes the data that are supplied are insufficient to identify a unique 
vibrating system; there is then a family of systems having the specified properties 
- an isospectral family. Chapter 7 describes how one can form such isaspectral 
families , and be sure that each member of the family has the necessary positivity 
properties. There are essentially two ways of forming families: algebraic, and 
differential. The former uses a carefully chcsen rotation to go from one member 
to another. The latter uses the idea of isospectral flow; a matrix can flow, under 
axalled Toda flow along a path ao that it retains the same eigenvalues and at 
the same time retains a particular structure and partkular positivity properties. 

Chapter 8 is concerned writh one particular type of vibrating system: a beam 
vibrating in flexure. This problem had been a severe stumbling block in the 
early history of inverse problems. 

Chapter 9 completes the first part of the bio ok with a study of modes, i.e. t 
normal modes, and nodes. This analysis depends heavily on the positivity study 
of Chapter 6. 

The second part of the book, Chapters 10-14, is concerned with continuous 
systems. The problems appear in two related forms, differential equations and 
integral equations. The integral equations, w r hich use the Greens function for 
the system, are the easier to analyse, for it is the Green’s function, Gantmacher 
and Krein’s kernel , that has the all-important positivity properties. Moreover, 
the Green’s function operator appearing in the integral equation is a concrete 
example of a positive compact self-adjoint operator in a Hilbert space, so that 
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we may immediately make use of t be well-developed theory of such operators, 
as described in Chapter 10. 

Chapter 11 uses this theory, and the fundamental Gel’fand- Levitan transfor- 
mation operator, to provide solutions to some inverse problems for the Sturm- 
Liouville equation. This equation, w'hich appears in three related forms, is the 
governing equation for the vibrating string and rod. The Chapter describes 
the classical approach, as well as some recent techniques that are more readily 
adaptable to computation. 

Chapter 12 discusses families of isxxspectral continuous systems. Chapter 
13 applies the Gel’fand-Levitan transformation to the inverse problem for the 
continuous Euler-Bernoulli beam. 

Chapter 14 is a short (too short) study of inverse nodal problems. While 
it is difficult in practice to measure a vibration mode, it is comparatively easy 
to locate the nodes of a particular mode. There is now’ a considerable body 
of research, due primarily to McLaughlin and Hald, that focuses on what nodal 
data is sufficient to identify, say, the mass distribution on a vibrating string, 
rod, or membrane, and how one can construct such a vibrating system horn a 
knowledge of some nodes of some modes. Section 14.4 briefly reports on this 
research. 

The book concludes with another short chapter on damage identification. 

The history of mathematics and the physical sciences leads to an important 
far-reaching conclusion: the study of one topic can throw light on many other 
topics, even on some which at first seem have no connection with the original 
topic. The study of inverse problems in vibration provides a clear example of 
this connectedness. On the one hand, there are topics in inverse problems that 
are illumined by knowledge in other fields, notably linear algebra and opera- 
tor theory; on the other hand the study of inverse vibration problems throws 
light on the classical direct problems by highlighting the fundamental qualitative 
properties of solutions. 

A remark on the quotations from Pascal’s Pens&es is in order. I used the 
translation by W.F. TVotter that appeared in Everyman’s Library, published by 
J.M. Dent & Sons in 1956. My copy is dated 2Gth April 1957 and contains an 
8d (old pence) ticket for the London Transport bus No. 73 from Euston Road to 
Stoke Newington, reminding me that the Pens&s were my daily bus reading to 
and from my 'digs’ when I was Assistant Lecturer in Mathematics at University 
College London. I chose the Penstes for the chapter captions because it is clear 
from his writings that Pascal considered the search for (Jod to be an inverse 
problem. His comments on the place of reason, heart and will in seeking a 
solution of the problem, though sometimes enigmatic, are as deep and relevant 
in 2004 as they were in 1651. I hope that these excerpts from the Pens&es will 
whet readers’ appetites for Pascals writings. 

The caption for Chapter 11 reminds me that many people have contributed 
to this book. Some were acknowledged in the Preface to the first edition. This 
new edition contains material taken from papers written with graduate students 
Brad Willrns, Mohamed Movahheddy, Hongmei Zhu and with colleagues Brian 
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Davies, Josef Leydold, Peter Stadler and Antonino Morrassi. In addition to 
these. I have freely taken from papers by numerous colleagues worldwide, as 
referenced in the bibliography. 

Parts of the book were read at the proof stage by Antonino Morrassi, Maeve 
McCarthy. Oscar Rojo and Michele Dilena. I thank them for pointing out many 
errors and shortcomings, some of which I have managed to correct. 

The book was typed by 'Racy Taves. Thank you for your stamina and your 
attention to detail. Colin Campbell helped us out with his understanding of 
the idiosyncracies of LaTeX. 

Finally, I acknowledge the patience and understanding of my wife, Joyce, 
who saw me immersed in books in my study for years on end. 

George Carrie* once remarked that the aim of mathematics is insight, 
numbers. It is the author’s wish that this book will provide insight into 
many interconnected topics in mathematics, physics and engineering that 
pear in the study of inverse problems in vibration. 



GALL. Gladw^U 
Waterloo. Ontario 
March, 2001 
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Chapter 1 

Matrix Analysis 



It is a bad sign when, on seeing a person, you remember his book. 1 

Pascal’s Pens^es 



1 . 1 Introduction 

The book relies heavily on matrix analysis. In this Chapter we shall present 
the basic definitions and properties of matrices, and provide proofs of some 
important theorems that will be used later. Since matrix analysis now has an 
established position in Engineering and Science, it will be assumed that the 
reader has had some exposure to it; the presentation in the early stages will 
therefore be brief. The reader may supplement the treatment here with standard 
texts. 



1.2 Basic definitions and notation 

We use the word iff to mean ‘if and only if’. A matrix is a rectangular array 
of real or complex numbers together with a set of rules that specify how the 
numbers are to be manipulated. 

A matrix A is said to haw order m x n if it has m rows and n columns. 
The set of all real matrices, i.e. t matrices with real entries, of order mxn, is 
sometimes denoted by R m x n . Following Horn and Johnson (1985) (183], we 
use the simpler notation M m , ni and say A € A / m .«. We write 



1 Blanc Pascal (1623-1662) lired among tbc French intelligcntia. and in that context it 
a bod sign; one should be known for more than just a book one hod written. When the first 
edition of this book was being translated into Chinese, tbc translator objected, for in 20th 
century China, it would be a goerf sign. If you met someone you knew who had written a 
book, you would mention it immediately! 
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The entry in row i and column j is a tj , and A is often written simply as 

A = (atf). 

Two matrices A, B are said to be equal if they have the same order m x n, 
and if 

= &«. (« = j = 1,2,..., n); 

Then we write 

A = B 

The transpose of the matrix A is the n x m matrix A r , whose rows are the 
columns of A. We note that the transpose of A r is A; we say that A and A r 
are transposes (of each other), and write this 

(A*)' = A. 

For example 



A = 


II 

1 

M <C 
-H CC 


1 2 ' 
2 G 




L 


—4 7 



are transposes. 

If m = n then the m x n matrix A is said to be a square matrix of order n: 
A € A/,,.,,; we abbreviate A/ n<r to A/„; thus A € M„. A square matrix that is 
equal to its transpose is said to be symmetric ; in this case 

A = A r , 

or alternatively 

= <*>». (i,j = 1,2, 

The set of real symmetric matrices of order n is denoted by S„. The matrix 

[ 1 2 9 ' 

A = I 2 •! G 

[ 9 G 3 

is symmetric. The square matrix A is said to be diagonal if it has non-zero 

entries only on the principal diagonal running from top left to bottom right. We 
write 

= diag(a n ,a 32 ,..,a m )- 



A = 



“li 

0 

0 

0 



0 

0 



0 

0 



0 

0 

0 



A = 



an oi2 
021 022 



Oml O r «2 



I 
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The imtf matrix of order n is 

I = I„= Av(l,l,...,l). 

The elements of this matrix are denoted by the Kronecker della 

The zero matrix of order mxnis the matrix with all its m x n entries zero. 

A matrix with 1 column and n rows is called a column vector of order n, and 
is written 

*1 
*2 

*« . 

The set of all such real vectors constitutes a linear vector space that we denote 
by Vn. 

The transpose of a column vector is a row vector, written 

* T = (*!»*■. •••.*■!■ 

Two matrices A, B may be added or subtracted iff they haw the same order 
m x n. Their sum and difference are matrices C and D respectively of the same 
order m x n, the elements of which are 

Cij = a,, + b,j , d,j = a,, - b,j. 

We write. 

C = A + B , D = A - B. 

The product of a matrix A by a number (or scalar ) k is the matrix kA with 
elements far,,-. 

Two matrices A and B can be multiplied in the sense AB only if the number 
of columns of A is equal to the number of rows of B. Thus if A has order mxn, 
B has order n x p then 

AB = C, 

where C has order m x p. We write 

A(m x n) x B(n x p) = C(m x p). (1.2.2) 

The element in row i and column j of C is c,j , and is equal to the sum of the 
elements of row i of A multiplied by the corresponding elements of column j of 
B. Thus 

N 

Cij = Oiibij + 4* . . . + Oinkij = ^2 Qik hii (1.2.3) 
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and for example 



r 2 3 1 1 


i 


2 


1 


0 ' 


[ 17 11 


l« • ’J 


0 

-1 


1 

0 


-1 

2 


0 

1 

■ 


- [ -6 8 9 7 



The most important consequence of this definition is that matrix multiplication 
is (in general) non-commutative, i.e., 

AB ^ BA. 

Indeed, if A is (m x n) and B is (n x p) then BA cannot be formed at all unless 
m = p. Even when m = p, the two matrices are not necessarily equal, as is 
shown by the example 



(1.2.4) 



In addition, this definition implies that there are divisors of zero; i.e., there can 
be non-zero matrices A, B such that 



A = 


1 

0 


-!]• 


B = 


1 

-1 


AB = 


[? 


jl 


BA = 


‘ 1 




1 1 


J 




-1 



AB =0. 



An example is provided by 



[ 1 1 1 1 1 1 _ r 0 °1 

[ 2 2 -1 -1 J - [ 0 0 J • 



The product of A (m x n) and a column vector x(n x 1) ts a column vector 
y(m x 1 } with elements 



3A = o.iZi+ 0 * 2*2 4*... + Oi„2:„ , (i = 1,2 
This means that the set of m equations 

011 *! + 012*2 + • • • + 01 r*n = 01 . 

021*1 + 022*2 + • • • + a 2n x n = y 2 . 



(1.2.5) 



( 1 . 2 . 6 ) 



a m ,x, + a m2 z 2 + . . . + a m „x n = n*,, 
may be written as the single matrix equation 



“li 

Oj, 



a, 2 

0,2 



“ml “m2 



Oln 

O-lr. 

“m~ 

Ax = y. 





» • 
^1 




yi 




x 2 








• 








. *» . 




. y<» . 



(1.2.7) 



( 1 . 2 . 8 ) 
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The product of an (n x 1) column vector x and its transpc^e x r (l x n) is an 
n x n symmetric matrix 



xx T = 



x\ X\X2 
X2X\ 



*1 



X\X n 
X2 X n 



(1.2.9) 



lull Z„X 1 - • - X 2 n 

On the other hand, the product of x T (l x n) and x(n x 1) is a (1 x 1) matrix, 
i.e., a scalar 

x t x = xJ + j3 + ... + j£. (1.2.10) 

This quantity, which is positive iff the x, (assumed to be real) are not all zero, 
is called the square of the Lj norm of x, i.e., 

INI 2 = x T x , ||x|| = (ij + x\ + . . . +a£)K 

The scalar (or dot) product of x and y is defined to be 

x T y = y T x =!,!(, + x 2 y 2 + . . . + J„y„. 

Two vectors are said to be orthogonal if 



x r y = 0. 



( 1 . 2 . 11 ) 

( 1 . 2 . 12 ) 

(1.2.13) 



It has been noted that matrix multiplication is rum- commutative. This holds 
even if the matrices are square (see (1.2.4)) or symmetric, as illustrated by 

r 1 2 1 r 1 - 1 1 = r 1 1 r 1 -‘if 1 2 l = f - 1 0 

[2 2 J [ — 1 lj [ 0 0j’[-l 1 J [ 2 2 J [10* 

(1.2.14) 

This example, which shows that the product of two symmetric matrices is not 
(necessarily) symmetric, hints also that there might be a relation between the 
products AB and BA. This result is sufficiently important to be called: 



Theorem 1.2.1 

(AB) T = B r A T , 

so that ui hen A, B. ore symmetric, then 

(AB) T = BA. 



(1.2.15) 



(1.2.16) 



Proof. Consider the element in row t, column j on each side of (1.2.15). 
Suppose A is (m x n), B is (n x p)> then AB ismxp and (AB) T is p x m. 
Then 



((AB) T )y = (AB)* = £ a )k b k „ 



k - 1 



and 



(B T A T )ij = ( row i of B r ) x ( column j of A r ) 
= ( column t of B) x ( row j of A) 

= EZ-, **«* ■ 
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1. If 



A = 



1 

2 

3 



find a square matrix B such that AB = 0. Show that if a 33 is changed 
then the only possible matrix B would be the zero matrix. 

2. Show that, whatever the matrix A, the two matrices AA r and A r A are 
symmetric. Are these two matrices equal? 



3. Show that if A, B are square and of order n, and A is symmetric, then 
BAB r and B r AB are symmetric. 

4. Show that if A. B. C can be multiplied in the order ABC. then ( ABC1 T = 
C r B T A T . 



5. If x is complex, then its L 2 norm is defined by 

INI 2 = |*i I 2 + M 2 + • - • + M 2 . 

Show that 

||x|| 2 = x*x 

where x* = Jt T , the complex conjugate transpose of x. 



1.3 Matrix inversion and determinants 



In this section we shall be concerned almost exclusively with square matrices. 
The determinant of a (square) matrix .4, denoted by det(A) or |A|, is defined 
to be 

det(A) = |A| = ± a Ul a 3ii • • • a„,„ ; (1.3.1) 

where the suffices are a permutation of the numbers 1,2,3, 

the sign is + if the permutation is even, and — if it is odd, and the summation is 
carried out over all n! permutations of 1,2,3, ... ,n. We note that each product 
in the sum contains just one element horn each row and just one element from 
each column of A. Thus for 2 x 2 and 3x3 matrices respectively 



a n o 12 
021 022 

Oil 012 
021 022 
031 032 



= OnO^ 




-o 12 o 2 i, 

011022033 

—011023032 



(1.3.2) 

Ul2<*23<*31 + <U3<*21<»32 

ai 2 U 21 du — H13U22G31 



The permutation *i,*2. is even or odd according to whether it may be 
obtained from 1 , 2 , . . . , n by means of an even or an odd number of interchanges. 
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respectively. Thus 1,3, 2, 1 and 2, 3, 1,4 are respectively odd and even permuta- 
tions of 1,2, 3, 4 because 

(1.2.3.4) — (1,3, 2, 4), 

(1.2. 3. 4) -(2, 1,3, 4) -(2, 3, 1,4). 

We now list some of the properties of determinants. 

Lemma 1.3.1 // too rotes (or columns) oj A ore interchanged, the determinant 
retains its numerical value, but changes sign. 

If the new matrix is called B then 

f>u = = a u , bj, = ajt , (j = 3, 4, . . . ,n) 

and 

det(B) = £ ± friufaafai '-fcnt., 

= 51 ± <22i a ai,- a a3i s * • * a«u ? 

= 51 03** - • • Oni. . 

But if *i>t2«i3>-*- s in is even (odd) then *2»*it*3>* • *i*n is odd (even), so that 
each term in det(B) appears in det(A) (and vice veisa) with the opposite sign, 
ao that det(B) = - det(A). 

Lemma 1.3.2 If (wo rows (columns) of A arc identical then det(A) = 0. 

If the two rows (columns) are interchanged, then, on the one hand, det(A) 
is unchanged, while on the other. Lemma 1.3.1, det(A) changes sign. Thus 
det(A) = -det(A) and hence det(A) = 0. 

Lemma 1.3.3 If one row ( column ) of A is multiplied by k then the determinant 
is multiplied bp k. 

Each term in the expansion is multiplied by k. 

Lemma 1.3.4 If two rows (columns) of A are proportional, then det(A) = 0. 
This follows from Lemmas 1.3.1, 1.3.3. 

Lemma 1.3.5 If one row (column) of A is added to another row (column) then 
the determinant is unchanged. 

If the matrix B is obtained, say, by adding row 2 to row 1 then 
bu = oii +02, * bji = <h* > i = 2,3 

Thus 

det(B) = £±*144**14--'**. = 

= E ± (°U. + °2i, ) a *j a 3., • • • ".u.. 

= E± a i<i °2i* a 3*, ■ ■ ■ a nin ± £ a Ml <h>,as<, ■ ■ • , 

and the first sum is det(A) while the second, having its first and second rows 
equal is zero. 
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Lemma 1.3.6 If a linear combination of rows (columns) of A is added to an- 
other row (column) then the determinant is unchanged. 

This follows directly from Lemma 1.3.5. We may now prove 

Theorem 1.3.1 If the rows (columns) of A are linearly dependent Mendet(A) = 
0. 

Proof. Denote the rows by af , aj, .... aj*. By hypotheses, there are scalars 
C\ , c 2 , . . . , c^ not all zero , such that 

ciaf + C2tg + ■ - • + c n a£ = 0. 

There is a c,- not zero; let it be Then 

= 5Z<c*/C".) a T- 

ial 



If the sum on the right is added to row m of A. the new matrix has a zero row, 
ao that its determinant, which by Lemma 1.3.6 is det(A), is zero ■ 

Before proving the converse of this theorem, we need some more notation. A 
minor of order p of a matrix A is t he determinant of a (square) submatrix of A 
formed by taking elements from p rows t|, t 2 , . . . >i p and p columns ji, j 2> . . . , j p . 
We denote the minor by 



Thus if 




i,»a> 

2 13 ' 
A = I -1 2 -1 
1 0 7 



(1.3.3) 



A(l;l) = 2, A(l,2; 1,2) = 



2 1 

-1 2 



= 5, A(l,2;2,3) = -2. 



There is an important special case. The minor of order n— 1 obtained by deleting 
the tth row and jth column of A is denoted by dy. Thus for the A in (1.3.3), 



dn = 



21 m A _1 4 
0 7 =14 ’ ai2= 17 



= -11, dia = 



-1 2 

1 0 



= - 2 . 



The minors ^ occur in the expansion of a determinant: for the determinant in 
(1.3.2) we may write 



det(A) = on (022033 — 023032 ) — ** 12(021033 — 1 



— °iidll — 012012 + 013013 



l) + 013(021032 “ 02203l) 
(1.3.4) 




J. Matrix Analysis 



This is called the expansion of det(A) along the first row. Thus for A in (1.3.3) 
we have 

33 = 2 x 14 - 1 x (-11) + 3 x (-2). 

The coefficients An,— 612,613 in (1.3.1) are called the cofactors of an, < 112,413 
respectively, and are denoted by An, A J2 , A 13 respectively. Thus we write (1.3.4) 
as 

det(A) = anAii +ai2Ai2 4- ai3Ai3 
an ai 2 ai3 
021 022 023 

031 032 O33 

If w»e take the cofactors of the first row and multiply them by the elements of 
another row. say the second row, then we get zero: 

021 022 023 
021 022 023 

031 032 O33 

The determinant on the left is zero because it has two rows equal. These two 
results, (1.3.5) and (1.3.6), are special cases of 

Theorem 1.3.2 

n 

Y “i * A Jk = det (A )6ij 
1-1 
n 

yi AktQkj = det(A)<5tf 
1-1 

uv/iere Sij is defined in (1.2.1). 

Proof. When i = j, so that 6 a = 1, these equations merely state the 
definition of a cofactor. When i ^ j they state that the determinant of a matrix 
with two row’s (or columns) equal, is zero ■ 

Now* compare equation (1.3.7) with (1.2.3). If we define a matrix B such 
that 

b kj = A jk (1.3.9) 

then we can write (1.3.7) as 

n 

Y = det(A)^ (1.3.10) 

1-1 

which, in matrix terms, states that 

AB = det(A)I (1.3.11) 

Likewise, (1.3.8) may be written 



(1.3.7) 

(1.3.8) 



= <»2lAu + a22i4l2 + 023^13 = 0 (1.3.6) 





BA = det(A)I. 



(1.3.12) 
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The matrix B is called the adjoint (or adjugate) of A and is denoted by adj(A). 
Thus equation (1.3.11), (1.3.12) state that 

A adj(A) = ad; (A) A = det(A)I. (1.3.13) 

We are now in a position to prove the converse of Theorem 1.3.1, namely 

Theorem 1.3.3 // det(A) = 0, then the rouw (columns) oj A are linearly de- 
pendent. 

Proof. We prove the result for the columns. That for the rows may be 
proved likewise. We will prove it by induction on n. It certainly holds, trivially, 
when n = 1, for then det(A) = a tl . Let a,, a^ . . . , a* be the columns of A, and 
suppose det( A) = 0. Either each set of n — 1 vectors selected from aj , a 2 , . . . a n 
is a linearly dependent set, in which case the complete set is linearly dependent 
as required, or there is a set of n — 1 vectors, which without Ices of generality 
we may take to be a! , a 2 , . . . a,,_ j , which is linearly independent. Now imagine 
creating a set of vectors bj, b 2 ,. . . , b„_i by deleting the ith row of each of the 
vectors ai,a 2 ,...,a n _i. For at least one value of i the set bj,b 2 , ...,b„_i 
must be linearly independent. By the inductive hypothesis, the (n — 1) x (n — 1) 
determinant formed from these vectors must be non-zero; at least one of the 
terms b*, in equation (1.3.10) is non-zero. If det(A) = 0, equation (1.3.10) 
states that 

n 

= 0 i,j= l,2,...,n (1.3.14) 

*-l 

Since a* = {ai*,a2*, . .. we may write the n equations (1.3.14) obtained 

by taking j = 1 , 2 ,..., n, as 

n 

X>»,a*=0 j = 1,2 n. (1.3.15) 

*-l 

For at least one value of j, not all the are zero; the columns a l? a 2 , . 
are linearly dependent ■ 

Theorem 1.3.4 The matrix equations 

Ax = 0. y r A = 0 

have non- trivial solutions x amt y respecth'ely iff det( A) = 0. 

Proof. The theorem is a corollary of Theorem 1.3.3. If ai , a2, . . . . an are 
the columns of A. then 

Ax = (ai,a 2 an]{xi,xa,...,zn} 

= nai + x 2 a2 + h x a a*. 

We can find xi, . . . ,x„, not all zero, such that 



xiai + x 2 aj H + j„an = 0 
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ill a], a],..., a,, are linearly dependent. By Theorem 1.3.3 this happens iff 
det(A) = 0. This happens in turn iff the row's of A are linearly dependent, i.e., 
y r A = 0 has a non-trivial solution ■ 



Theorem 1 . 3.5 // A.B are square mafrtees of order n then 

det(AB) = det(A).det(B) 

The proof of this result is left to Ex. 1.3.3. 

The square matrix A is said to be singular if det(A) = 0, non-singular or 
invertible if det(A) ^ 0. Theorem 1.3.4 shows that if A is non-singular, then 
the equation Ax = 0 has only the trivial solution x = 0. Ex. 1.3.5 extends this 
result: if A is non-singular, then the matrix equations AS = 0, TA = 0 have 
only the trivial solutions S = 0, T = 0; when A is non-singular there are no 
divisors of zero. 

The matrix R is said to be an inverse of A if AR = I. 

Theorem 1.3.6 If A has an inverse, it is unique, and RA = I. 

Proof. Suppose AR = I. Theorem 1.3.5 shows that 

det(A).det(R) = det(I) = 1 (1.3.16) 

so that det(A) 0: A is non-singular. If Ri,Rj were two inverses, then 
AR, = 1 = AR 2 , so that AIR, - Rj) = 0. But A is non-singular, so that 
R, — R 2 = 0: Rj = R,. Now if AR = I then ARA = A, i.e., A(RA - 1) = 0. 
But A is non-singular, so that RA -1 = 0, i.e.. RA = I ■ 

Theorem 1.3.6 shows that if A has an inverse, then A is non-singular. The 
logical negative of this statement is that if A is singular it does not have an 
inverse. We now prove the converse. 

Theorem 1 . 3.7 If A is non-singular, then it has an inverse. 

Proof. If A is non-singular, then det(A) ^ 0, and equation (1.3.13) may be 
written 

AR = RA = I, (1.3.17) 

where R= adj(A)/det(A) ■ 

If A is non-singular, its unique inverse is denoted by A -1 . We have 

AA 1 = A -I A = I. (1.3.18) 

Theorem 1 . 3.8 The equation 

Ax = b (1.3.19) 

either has a unique solution, if A is non-singular: or if A is singular, it has a 
solution only for certain b. 
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Proof. If A is non-singular then 

x = A" '(Ax | = A _, b 

is the unique solution. If A is singular, then there is one (or more) y such that 

y T A = 0 . 

Then 

y T (Ax) = y T b = 0 

so that (1.3.19) has a solution only if b is orthogonal to any y which satisfies 
y 1 A = 0 . If A is singular then Ax = 0 has one or more solutions Xi,X2,. . . ,x m , 
so that if xo is one solution satisfying Axo = b, then 

m 

x = x o + E C ‘ X > (1.3.20) 

<— 1 

is also a solution for arbitrary ci, ca, ... ,Cm ■ 

Note that trying to solve Ax = b by actually finding the inverse of A, 
is an extremely wasteful and clumsy procedure. Finding A " 1 is equivalent to 
living Ax = b for oil possible b, not just for the given b. Techniques for 
living Ax = b form the subject matter of numerical linear algebra, for which 
see Bishop, Gladwell and Michaelson (1965) (33) or Golub and Van Loan (1983) 
(135). Note also that we haw not in fact shown how* to find one solution Xq if B 
is in fact orthogonal to all solutions of y r A = 0 ; this too is covered in numerical 
linear algebra. 

In numerical linear algebra the starting point of almost ail the procedures 
for solving linear equations such as (1.3.19), whether A is square or not, or of 
finding determinants, is Gaussian elimination. This is a systematic reduction of 
an array (o^) to (usually) upper triangular form by subtracting multiples of one 
equation from another. Lemma 1.3.6 show's that the determinant of coefficients 
is unchanged under such an operation. 

The application of Gaussian elimination to the equations 

xi 4- 3x2 4- 2x3 = 6 

2 xi + 5x2 4- 6 x 3 = 13 
3xi 4- 4x2 4- 7x3 = 14 



would proceed as follows; only the coefficients need be retained: 



: 


3 


2 


6 


— * 


1 


3 


2 


6 


— ♦ 


1 


3 


2 : 


G 


2 


5 


G 


13 




0 


-1 


2 


1 




0 


-1 


2 : 


1 


3 


4 


7 


14 




0 


-5 


1 


-4 




0 


0 


-9 : 


-9 



The determinant of A is 1 x (-1) x (-9) = 9. The last of the new equations 
gives —9x3 = —9, X 3 = 1; when substituted in the new' second equation this 
gives -X 2 = 1 - 2 x 3 = -1, X 2 = 1; then xi 4-3 4-2 =6 gives xi = 1. 
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Exercises 1.3 

1. Show that if A is upper (lower) triangular, i.e., a,j = 0 if j > i(j < *), 
then 

det(A) = 011022... a m - 

2. If 

1 3 2 

2 5 6 

3 4 7 

find A" 1 . Verify that A A 1 = A 1 A = I 

3. Prove that if A, B are square matrices of order n, then 

det(AB) = det(A).det(B). 

Hint: consider the 2n x 2n matrix 




Show that det(C) = det(A). det(B). Now subtract multiples of rows (n+1) 
to 2n from rows 1 to n to delete all elements in the top left quarter of C. 
The elements in the top right quarter will be those of — AB. 

4. Use Gaussian elimination to solve the equations 

xi + 2x 2 + 4x3 + 8x* = -9 

X2 + 3X3 + 2x« = 1 

xi + 2X2 + 5x3 + Gx« = -3 

-XI + 3X2 + <1X3 + 7x« = -10 

5. Show that if A is non-singular, then the matrix equations AS = 0 and 
TA = 0 have only the tririal solution S = 0, T = 0. respectively. 

1.4 Eigenvalues and eigenvectors 

If A and C are square matrices of order n then the equation 

Cx = A Ax (1.4.1) 

will have a non-trivial solution x (i.e., one for which ||x|| ^ 0) iff the matrix 
C - AA is singular, i.e., the scalar A satisfies the determinantal or characteristic 
equation 

det(C - AA) = 0. (1.4.2) 

The roots of this equation are called the eigenvalues of the matrix pair (C T A); 
they may be real or complex. If A is an eigenvalue, a vector x satisfying (1.4.1) 
is called an eigenvector corresponding to A. 
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In many mathematical texts, attention is focused almost exclusively on the 
case when A = I. In this case A is said to be an eigenvalue of C. The problem 
(1.4.1) is called the generalised eigenvalue problem. In Mechanics there are 
many problems in which two matrices, C, A appear, and it will be convenient 
to develop the theory for this case. 

The eigenvalue theory for general, i.e., not necessarily symmetric matrices 
C, A, is extremely complicated. (See Ex. 1.4.8). However, for all, or almost all, 
the problems encountered in this book, the matrices C, A have special properties: 
they are real and symmetric, and one at least is positive definite , defined as 
follows. 

Suppose A is real and symmetric, and x is a real nxl column vector. The 
quantity x r Ax is a scalar. Written in full it is 

x r Ax = o„x? +2a u x l z a +- • ■+2a ln x l x n +a„x\ + • • ■+2a 2n x 3 x n +• • +a nf ,x^. 

(1.4.3) 

This is called a quadratic form. In many physical applications the kinetic energy 
and the potential energy of a mechanical system may be expressed as quadratic 
forms in the generalised velocities or displacements, respectively. The kinetic 
energy of a system is always positive, unless all the generalised velocities are 
zero. This leads us to a definition. The matrix A is said to be positive definite if 
| |x|| ^ 0 implies x T Ax > 0. (Clearly, if ||x|| = 0, so that x\ = 0 = f2 = . . 
then x T Ax = 0.) If A satisfies the weaker condition, that ||x|| ^ 0 implies 
x r Ax > 0, i.e., there is a vector x such that ||x|| / 0 and x T Ax = 0, then A is 
said to be positive semi-definite. We will find later that the matrix associated 
with the potential energy of an unanchored system is positive semi-definite; there 
is a vector x corresponding to a rigid body displacement of the system, for which 
the potential (or sfrotn) energy is zero. 

Theorem 1.4.1 7/C. A are real and symmetric, and A is positive definite . then 
the eigenvalues and eigenvectors of (1.4 -1) are real. 

Proof. Suppose A. x possibly complex, and with ||x|| ^ 0. are an eigenpair 
of (1.4.1), multiply both sides by x* = X r to obtain 

x # Cx = Ax* Ax (1.4.4) 

The quantities x* Ax and x*Cx are both real. This is ao because x*Ax, for 
instance, is a scalar, and therefore equal to its ow r n transpose. Thus 

a = x*Ax = (x*Ax) r = x T A r x = x r Ax = (x*Ax) = a 

but if a = a, then a is real. Similarly, x*Cx ts real. Moreover, if ||x|| 0, i.e., 

at least one element in x is not zero, then a is strictly positive, i.e., a > 0. For 
let x = u -f tv where u, v are real, then 



x* Ax = (u T - »v r ) A(u + tv) = u r Au 4- t{u r Av - v r Au) + v T Av. 
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But since x*Ax is real, the imaginary term is zero, and thus 

x* Ax - u T Au 4- v T Av > 0. 

The inequality is strict because either at least one element of u is non-zero, in 
which case u T Au > 0; or if u = 0. at least one element of v is non-zero, in 
which case v T Av > 0. 

Now return to equation (1.4.4); x*Cx and x*Ax are both real and x*Ax is 
positive. Hence 

A = x # Cx/x # Ax 

is real. Since A is real, the vector x. obtained by solving a set of simultaneous 
linear equations with real coefficients, is real. Therefore, x # = x T , and we can 
WTite 

A = x r Cx/x r Ax. ■ 

This ratio is often called, and we will call it. the Royleigh Quotient corre- 
sponding to equation (1.4.1). (It was Lord Rayleigh (Rayleigh (1S94) (290)) 
who, in his classical treatise Theory of Sound used this quotient to take the first 
steps towards the variational treatment of eigenvalues. We discuss this further 
in Chapter 2.) We write 

Ah = R(x) = x T Cx/x T Ax (1.4.5) 



Ex. 1.4.7 shows the necessity of having one of the matrices A. C. positive 
definite. 

The conditions which must be satisfied if a (symmetric) matrix A is to be 
positive definite or positive semi-definite may be expressed in terms of the prin- 
cipal minora of A. A principal minor of order p of a matrix A (symmetric or 
not) is a determinant of a submatrix formed from p rows ii.ij, . . . ,i P and the 
same p columns Thus for A in (1.3.3), 



2 1 

-1 2 



2 3 2 1 

1 7 ' 0 7 



2 1 3 
-1 2 4 
1 0 7 



are all principal minors. In the notation of Section 1.3, a principal minor is 

There is a special notation for the leading principal minora of A, these are 
as follows: 



D, = a n , D 2 = 



an 

aji 



D n = |A|= det(A). (1.4.6) 



Now we may state 



Theorem 1.4.2 The symmetric matrix A is posit tie definite iff the leading 
principal minors (D,)i ait alt positive. A icill be positive semi-definite iff 
(Di)l~ l >0, D„= 0. 
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This will not be proved until Chapter 5. Note that since D„ = det(A), this 
states that a positive-definite matrix is non-singular, and a positive semi-definite 
matrix is singular. 

We may now refine Theorem 1.4.1 to give 

Theorem 1.4.3 If C, A ore real and symmetric and A is positive definite then 
equation (1.4.1) will hate n real eigenvalues, although they need not be distinct. 
If C is positive definite they will be positive , if C is positive semi-definite they 
will be non-negative. 

Proof. Equation (1.4.2) may be expanded in terms of the coefficients Cy — 
X Oy \ the result is an nth degree polynomial equation for A, namely 

A(A) = det(C - AA) = A 0 + A, A + A*A 2 + • - • + A„A n = 0. (1.4.7) 

Most of the coefficients A, are complicated functions of a i j and Cy, but the first 
and last may be easily identified, namely 

A 0 = det(C), A„ = (-1)" det(A). (1,1.8) 

Since A is positive-definite. det(A) > 0 so that A,, ^ 0. This means that 
equation (1.4.7) is a proper equation of degree n with n roots (At)". This proves 
the first part of the Theorem. If C is positive-definite, then both numerator and 
denominator of the Rayleigh Quotient (1,1.5) will be positive, so that (At)? > 0. 
If C is only positive semi-definite, then the numerator of the Rayleigh Quotient 
is only positive or zero (i.e., non-negative), so that the A, are non-negative. 
Moreover, since AiA2...A n = (-)"Ao/A n = det(C)/det(A) equation (1,1.7) 
will have at least one zero root when det(C) = 0 ■ 

Under the conditions of Theorem 1.4.3 the eigenvalues (A,)? may be labelled 
in increasing order: 

0<Ai ^A 2 <...<A n . (1,1.9) 

Theorem 1.4.4 Eigenvectors Ui.u j corresponding to two different eigenvalues 
£ Xj) of the symmetric matrix pair (C, A) are orthogonal w.r.t. both 
A and to C, i.e,, 



ufAuj = 0 = u[ Cuj. ( 1 , 1 . 10 ) 

Proof. By definition 

Cu, = AtAu., Cu j = AjAiij. (1.4.11) 

Transpose the first equation and multiply it on the right by uj^ ; multiply the 
second equation on the left by uj\ to obtain 
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Subtract these two equations to yield 

(A 4 -A > K T Au i =0. 

But Xi- Xj ^0 y so that uf Auj = 0, and hence uf Cu^ = 0. ■ 

Premult iplying equation (1.4.11) by we find 

Ci = u[Cu % = A.u^Au, = A %Qi (1.4.12) 

Sometimes, we will normalise an eigenvector m w.r.t. A; then <n = l,cj = A*. 
An important corollary of this result is 

Theorem 1.4.5 If the symmetric matrix paxr (C, A) has distinct eigenvalues 
(A«)| , and A is positive- definite, then the eigenvectors u* are tmeariy indepen- 
dent , and therefore span V Mf the space of n -vectors. 

Proof. The eigenvectors u* are linearly independent; for suppose 



multiplying by uf A we have 

Q,(u, r Au,) + Q 2 (u, r Auj) + • - • + OnluT Au„) = 0. 

But Ilf Au; = 0 if t ^ j, so that only the ith term in this equation is non-zero, 
and hence 

o,-(uf Au,) = 0. 

Since A is positive definite, uf Au, > 0 and a, = 0; all the (o,)" are zero; the 
u, are linearly independent. Any vector ueln may be written uniquely as 



n 

U = y^Q;U; (1.4.13) 

1-1 



where 

O; = uJ’Au/ujAU;. ■ (1.4.14) 

In this book we are not concerned with methods for computing eigenvalues 
and eigenvectors. A simple treatment of the classical techniques may be found 
in Bishop, Glad we 11 and Michaebon (19C5) [33). A comprehensive account of 
modern techniques is given by Golub and Van Loan (1983) (135). The classical 
treatise on the symmetric eigenvalue problem is Parlett (1980) [264). We are 
concerned only with the qualitative properties of eigenvalues. 



Exercises 1.4 
1. If 



A = 



1 -1 

-1 2 -1 
-1 2 -1 
-1 1 

show that A is positive semi-definite. For what x is Ax = 0? 
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2. Show that A -1 is positive definite iff A is positive definite. 

3. Verify the conditions given in Theorem 1 .4.2 for A to be positive definite, 
when n = 2, by writing 

x T Ax = a n xj + 2 a l2 iiX 2 + <h 2*2 

Extend the analysis to n = 3. 

4. Find the eigenvalues and eigenvectors of the pair 



2 


-1 


0 ' 




' 1 


-1 


2 


-1 


, A = 


1 


0 


-1 


2 




• 



Hint: replace the eigen-equation by the equivalent recurrence relation 
— i r _i + (2— X)x r — 2V+ 1 = 0 with appropriate end conditions for r = 1 , r = 
3, and seek a solution of the form x r = A cos rti+ B sin rO. Generalise this 
result. 

5. Show that if n = 3, A is symmetric, and D ly D 2 , D 3 of equation (1.4.6) are 
all positive, then all the principal minors of A are positive. Hint: write 
an det(A) as a 2 x 2 determinant with elements which are minors of A of 
order 2. This is a particular case of a general result, 9ee e.g., Gantmacher 
(1959) [971. 

6. Show that the real symmetric matrix A has positive eigenvalues iff it is 
positive-definite. 




The eigenvalues are not real. Where does the argument used in the proof 
of Theorem 1.4.1 break down? 



8. Take 




Show that equation (1.4.1) has only one eigenvalue and one eigenvector, 
so that the eigenvectors do not span the space Vz. This is the kind of 
difficulty attending the non-symmetric eigenvalue problem. 




Chapter 2 

Vibrations of Discrete 
Systems 



Our nature consists in motion; complete rest is death. 

Pascal’s Penates, 129 



2.1 Introduction 

The formulation and solution of the equations governing the motion of a discrete 
vibrating system, i.e., one which has a finite number of degrees of freedom, have 
been fully considered elsewhere. See for example. Bishop and Johnson (I960) 
(34), Bishop, Glad well and Michaelson (1965) (33), Metrovich (1975) (234]. In 
this chapter we shall give a brief account of those parts of the theory that will 
be needed for the solution of inverse problems. 

Throughout this book we shall be concerned with the infinitesimal vibration 
of a conaert'ative system about some datum configuration, which will usually be 
an equilibrium position. 

Before embarking on a general discussion we shall firet formulate the equa- 
tions of motion for some simple vibrating systems. 

2.2 Vibration of some simple systems 

Figure 2.2.1 shows a vibrating system consisting of n masses connected by lin- 
ear springs of stiffnesses (fc r )y . The whole lies in a straight line on a smooth 
horizontal table and is excited by forces (F r (t))y. 

Newton’s equations of motion for the system are 

tn.ilr = F r + 0,*, - 9 r , r = 1,2,. ..,n — 1, (2.2.1) 

nu.u,, = F„-e n , ( 2 . 2 . 2 ) 
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where * denotes differentiation with respect to time. Hooke’s law states that the 
spring forces are given by 

0 r = Mtir -tir-i) , r = 1,2,... t n. (2.2.3) 

If the left hand end is pinned then 

no = 0. (2.2.4) 

Forced vibration analysis concerns the solution of these equations for given 
forcing functions F r (t). Free vibration analysis consists in finding solutions to the 
equations which require no external excitation, i.e., F r (t) = 0, r = 1,2, . . . ,n, 
and which satisfy the stated end conditions. 




The system shown in Figure 2.2.1 has considerable engineering importance. 
It is the simplest possible discrete model for a rod vibrating in longitudinal mo- 
tion. Here the masses and stiffnesses are obtained by lumping the continuously 
distributed mass and stiffness of the rod. Equations (2.2.1) - (2.2.4) a lsx> describe 
the toreional vibrations of the system shown in Figure 2.2.2.. provided that the 
Uy. k T1 nv are interpreted as torsional rotations, torsional stiffnesses and mo- 
ments of inertia respectively. Such a discrete system provides a simple model 
for the torsional vibrations of a rod with a continuous distribution of inertia and 
stiffness. 



There is a third system which is mathematically equivalent to equations 
(2.2.1) - (2.2.4). This is the transverse motion of the string shown in Figure 
2.2.3 which is pulled taut by a tension T and which is loaded by masses (mr)*. 
(But note that the string showir in Figure 2.2.3 has its right hand end filed, 
rather than free, as in Figures 2.2.1 and 2.2.2. In order to simulate a string 
with a free end, the last segment of the string must be attached to a massless 
ring that slides on a smooth vertical rod.) If in accordance with the assumption 
of infinitesimal vibration, the string departs very little from the straight line 
equilibrium position, then the equation governing the motion of mass m r may 
be derived by considering Figure 2.2.4. 
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Figure 2.2.2 - A torsionally vibrating system 



m 




Figure 2.2.3 - n musses cm a taut string 




u 



f 



Figure 2.2.4 ~ forces acting on the moss m*. 
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Newton’s equation of motion yields 

fUriir = F, + rsinQr+i -Tsinor , (2.2.5) 

= F.+ffr, i-0,, (2.2.6) 

where, for small deflections, we may take sina r = a„ 

Or = Tor = kr(Ur - Ur-l) , *r = T/t r . 

In order to express equations (2.2.1) - (2.2.3) in matrix form we use (2.2.3) 
to obtain 

«M*r = F r + k rtl tV*l-(*,.l+fcrK + *rUr-l ,">1.4. = F n - k,,U,,+ , 

which yields 




This equation may be written 



Mii + Ku = F (2.2.8) 

where the matrices M, K are called respectively the inertia (or moss) and the 
stiffness matrices of the system. Note that both M and K are symmetric ; this 
is a property shared by the matrices corresponding to any conservative system. 
We note also that both M, K are positive-definite. In this particular example 
the matrix M is diagonal while K is tridiagonal, i.e., its only non-zero elements 
are on the principal diagonal, and the two neighbouring diagonals, called the 
codiagonals. 

Equation (2.2.3) can afc*> be constructed by introducing 0 = |#i, 02 , . . . ,0 n ) 
and noting that 




(2.2.9) 



which will be written 



0 = KE r u . 
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where 



E = 



■ 1 -1 


0 ••• 


0 ' 




1 


1 


1 •• 


- 1 ■ 


0 1 


-1 ••• 


0 




0 


1 


1 •• 


- 1 


0 


1 


-1 


,E-' = 


0 


0 


0 •• 


- 1 


0 ... 


0 


1 

* 




0 


0 


0 •• 


- 1 



( 2 . 2 . 10 ) 



and k = dtag{fc!,fca, 

Using the matrix E, we may write equation (2.2.1) - (2.2.2) in the form 



Mu = -E0 + F . 



so that on using (2.2.9) we find 



Mil + EKE r u = F . 



( 2 . 2 . 11 ) 



and ~ 

K = EKE . (2.2.12) 

For free vibration analysis there are two important end conditions. The right 
hand end may be free, in which case there is no restriction on the (uj)", or it 
may be fired, in which case ti„ = 0. 



Exercises 2.2 

1. Verify that the stiffness matrix in equation (2.2.7) satisfies the conditions 
of Theorem 1 .4.2. Obtain a proof that applies to principal minors of any 
order t, such that 1 < t < n. 

2. Consider the multiple pendulum of Figure 2.2.5. 




Figure 2.2.5 - A compound pendulum mode up of n inextensible strings 
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Show that the kinetic and potential energies of the system for small oscilla- 
tions are given by 

2 T = m,|tf + m 2 y| + . . . + 

2V = o,4 + g2 l»^-K-. + o n , *-f- ,>a 

where o T = g m,. 



2.3 TYansverse vibration of a beam 



Figure 2.3.1 shows a simple discrete model for the transverse vibrat ion of a beam; 
it consists of n + 2 masses (mr)-i linked by massless rigid rods of lengths (£ r )o 
which are themselves connected by n rotational springs of stiffnesses (Av)J . The 
mass and stiffness of the beam, which are actually distributed along the length, 
have been lumped at n + 2 points. 



The discrete system is governed by a set of four first-order difference equa- 
tions. which may be deduced from Figure 2.3.2. 




Figure 2.3.1 - .4 discrete model of a vibrating beam 
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Figure 2.3.2 - The configuration around m. 



For small displace in aits, the rotations are 

B, = (xi T -Ur-i)ltr , r = 0,l,...,n. 

If the rth spring has rotational stiffness A^, then the moment r 9 needed to 
produce a relative rotation 0 r + x —0 r of the two rigid rods on either side of m r 
is 

T r = *r4l(*r+I — *r), T = 0, 1,. . . ,f» - 1. 

Equilibrium of the rod linking m r and m r +i yields the shearing forces 
❖r = (r r - r r4l )/^ r+1 , r = — l,0,...,n - 1, 
while Newton's equation of motion for mass m, is 

mriir =4>r-4>r-i , r = -1,0,..., n. 

Here <^_ 2 , 4>„ r -i> T » denote external shearing forces and bending moments, 
respectively, applied to the ends. 

Suppose that the left hand end is clamped so that 

u-i=0 = uo, 

(m r )i move, and the governing equations may be written 
0 = L -l E r u, (2.3.2) 

r = KE r 0, (2.3.3) 



then only the 
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<*> = L-'Er - <^ I r B e B1 (2.3.4) 

Mu * —E# + ^ n e„, (2.3.5) 

where u= {«i,Ua,-.. ,11,], 9 = {9 ly 0 2 ,... .*«}. T = 

•>= K = <fiay(* r ), L = <ftay(f r ), M = <ftay(m r ), 

e„ = {0,0,.. .,0,1} and E is given in equation (2.2.10). 

Equations (2.3.2) - (2.3.5) may be combined to give 

Mii*Ku = ^e n + C , r„Ee nI (2.3.6) 

where 

K = EL-'EfCE r L-'E T . (2.3.7) 

This equation has the same general form as equation (2.2.8). We note that 
M and K are again symmetric and positive-definite, M being diagonal, and K 
being pentadiogonal. 



2.4 Generalised coordinates and Lagrange’s equa- 
tions: the rod 

The idea that a discrete system is one composed of a finite number of masses 
connected by springs is unnecessarily restrictive. The general concept is that 
of a system whcee motion is specified by n generalised coordinates ($ r )j that 
are functions of time t alone. The systems considered in Sections 2.2, 2.3 are 
indeed discrete in this sense and the generalised coordinates corresponding to 
the system in Figure 2.2.1 are . However, the more general concept would 
also cover, for instance, a model of a non-uniform longitudinally vibrating rod 
constructed by using the finite element method (see for example. Zienkiewicz 
(1971) (343)), Strang and Fix (1973) [311). 

In such a model, shown in Figure 2.4.1, the rod is first divided into n + 1 ele- 
ments. In the rth element, shown in Figure 2.4.2., the longitudinal displacement 
p(x,f) is taken to have a simple linear form. 



y(x,f) = Vr(f)(»-0 + yr4l(^ . *r<*<* r +l. (2-4.1) 



where 

i = ( X -Xr)//r, 

runs from 0 at the left hand end of the element to 1 at the right. Equations 
(2.1.1) with r = 0,1 ,..., n express the displacement at every point of the rod 
in terms of the n + 2 generalised coordinates (yr)o* 1 - When the end conditions 
are imposed there will be. as before, only n coordinates (pr)i . 
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Figure 2.4-1 - *4 rod divided into elements 




Figure 2.4-2 - One element of the rod 



When the finite element method is used, it is not possible to set up the 
equations of motion by using Newton's equations of motion, for there is no actual 
‘mass' to which forces are applied. Instead we may use Lagrange's equations. 
For a conservative system with kinetic energy T and potential or strain energy 
V . which are functions of n coordinates Lagrange's equations state that 

(r=1 ' 2 ">• < 21 - 2) 

Here Q r is the generalised force corresponding to qr in the sense that the work 
done by external forces acting on the system when the system is displaced from 
a configuration specified by (qr) \ to one specified by (^r + is 

n 

6W* = ^ QrSqr. 
r-1 
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For the system shown in Figure 2.2.1 the kinetic and potential energies are 

7 = < 2A3 > 



r -1 



r -1 



and Qr = F r («). Thus 

_ = m r y r1 — = -Mpr,l - Vr) + *V-l(V, ~ 1/r-l)- 
and equation (2.4.3) yields (2.2.1). 

For the finite element model of Figure 2.4.1, the kinetic and potential energies 
of the system will be 

T =\ f Q SpWM**, 

where S(x) % p(x) % E{ jr) are the (possibly variable) cross-sectional area, density 
and Young’s modulus of the rod. On inserting the assumed form of y{x*t) given 
in (2.4.1) we find 

S{ Xr + trt )p(Xr + *rO(Vr(l - 0 + Vr+ltftrdl, (2.4.4) 



V = \Y1) ( 5(*r + tri)E{x, + trOlVr.l ~ ») PC'*- (2-4.5) 

C)n carrying out the integrations, perhaps numerically if S(x),p(x),E(x) are 
variable, we may write 



, «♦ 1 


«4l 

X] m r .v,v.> 

J-0 


(2-4-6) 


v=L n £ 

If the rod is fixed at both ends, then 


j-o 


(2.4.7) 


tA> = 0 


= Vm+1> 


(2.4.8) 



**> that all the sums in (2.4.6), (2.4.7) run from 1 to n. In this case 
OT ^ . OV 



or V' - ov V* . 

857 = 2-™”*. g* =!>•*. 
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and equation (2.1.2) yields the following equation for free vibration: 

n n 

£ m '*V» + ]T] Av * f, * = 0, < r = l,2,...,n). 

»-l J-l 

This equation may, as before, be condensed into the matrix equation 

My + Ky = 0 (2.4.9) 

We note that, for the rod with the kinetic and potential energies given by 
(2.4.6). (2.4.7), the matrices M.K are symmetric, tndiogoruil matrices with 
sign properties. They are tridiagonal because m r »,k r t are zero unless r = s 
or r = s ± 1. The sign properties may be deduced from (2.1.4), (2.4.5): the 
codiagonal elements m rt r+i,mr,r-i of M are positive, while A>, r +i, Av.r-i are 
negative. Thus 





01 bi 




ci -di 


M = 


b, a 2 


,K = 


-<f, C3 




■■ '• bn-l 




-dn-l 




6n-l On . 




~d»-\ C„ 



(2.4.10) 

These sign properties of M. K will later be shown to have important implications 
for the qualitative properties of a vibrating rod. 

On the basis of these examples we now pass to the general case. For a con- 
servative system with generalised coordinates {Qr)\ w r hich specify small displace- 
ments from a position of stable equilibrium, the kinetic and potential energies 
will have the form ^ 

(2.4.11) 

r-1 .-1 

(2.4.12) 

r-1 — 1 

where the matrices M = (m™) and K = (Av*) are symmetric, in that 

m Mr = mrj. ktr = fcj. 

The equations governing free vibration may be written 

Mq + Kq = 0. (2.4.13) 

We note that equations (2.4.11), (2.4.12) may be written 

T=i^ T Mq, V=lq T Kq. (2,1.14) 



It ks not possible for any arbitrarily chosen symmetric matrix M to be an 
inertia matrix, because the kinetic energy T is an essentially p<»itive quantity, 
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i.e., it is always positive except when each of the q r is zero, in which case it is 
zero. Thus M must be positive definite (see Section 1.4). 

The restrictions on the matrix K are slightly less severe since, although the 
strain energy w ill always be positive or zero, it will actually be zero if the system 
has a rigid-body displacement. Notice, for example, that the V of (2.4.5) will 
be zero if y is the rigid body displacement 

V0=Pl =-•* = lfc = Vn+l- 

This will be a possible displacement of the system in Figure 2.2.1 only if both 
ends are free. We conclude that if the system is not constrained so that one 
point is fixed, then K is positive semi- definite. 



Exercises 2.4 



1. Use equations (2.4.4), (2.4.5) to evaluate the mass and stiffness matrices 
for a uniform rod in longitudinal vibration subject to the end conditions 
(2.4.8). 

2. Use the form (2.4.5) of the strain energy of the rod to show that the stiffaess 
matrix K for a rod fixed at the left and free at the right has the form 



* l +*2 

-*2 



-*2 

k? + k 3 —A 



K = 



-*«-i 

— fcw-l kr 



2.5 Vibration of a membrane and an acoustic 
cavity 

Over the last three or four decades, computational vibration analysis has de- 
veloped to such an extent that it can analyse the vibration of almost anything: 
rods, beams, plates, trusses, steel and concrete buildings, bridges, aircraft, and 
so on. Inverse vibration analysis in the strict form we consider in this book can 
hope to encompass only comparatively simple structures: strings, rods, beams, 
membranes and acoustic cavities and, even now, inverse problems for membranes 
and cavities are still open; all we can do is find some qualitative properties of 
the vibration. The vibrations of a membrane and of an acoustic cavity are 
mathematically similar: both involve just one scalar quantity, the transverse 
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displacement u(x, y ), for the membrane under unit tension; the excess pressure 
p(z,y,z), for the acoustic cavity. Both are governed by a wave equation 



. - - . d 1 & 

Au ~ p W A ~o^ + W' 

for a membrane with ma& density p{x t y), and 

a 2 . a 2 . a 2 



Ap - p W' A -p T ^’F 



(2.5.1) 



(2.5.2) 



for the acoustic cavity. 

To set up the finite element model FEM of a membrane we consider the 



T~l! {/***"• 

/ =lj j (vu) 3 dxdy. 



(2.5.3) 

(2.5.4) 



The simplest FEM is based on triangulation. For an arbitrary triangular element 
A, ft, ft as shown in Figure 2.5.1, we take 



u{x y y) =o + bz + ey 



(2.5.5) 




Figure 2.5.1 - .4 triangular finite element 



If u takes the values ui,tt2> **3 at the vertices Pi. ft, ft respectively, then 

ti, =a + hi, +c&, i = 1,2,3. (2.5.6) 

We can solve these equations for a, 6,c and hence express T, V for one element, 
i.e. t T ei V e , as quadrat ic forms 

r« = iurM.u« (2.5.7) 
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V, = ^ufK,n« (2.5.8) 

with coefficients which are functions of the coordinate (Zitft),* = 1,2,3. We 
are not particularly interested in the magnitudes of the coefficients; w r e are more 
interested in their signs. 

First we investigate the elements of K*. Equation (2.5.8) give 

bA = u I (ir 2 -y 3 ) + u 2 (if 3 -y I ) + u 3 (i( 1 -jfc) 
cA = U,(x 2 -X 3 ) +ti2(23-I,) + U 3 (j 1 - Z 2 ) 



where 



A = 



1 H la 
1 i2 la 
1 zj K3 



= 2 area(ftftft). 



Since (yu) 2 = 6 2 + c 2 , the coefficient of, say, uit*2 in V, is 

-{(is - ii)(is - 12) + (Vs ~ Vi)(to ~ to)}/|A| = -|PiP 3 |.|P 2 Psl o«VIA|. 



Users of finite element methods have found that compact, i.e., acute angled, 
triangles give more accurate computational results than elongated triangles that 
haw an obtuse angle. 

If the triangle has all its angles acute, then A'u, and A - 31t are all neg- 
ative: K. has the sign pattern 



K, = 



+ 



+ 




(2.5.9) 



To find the signs of the coefficients in T e , it is convenient to WTite (2.5.7) 
in terms of the areal coordinates , y), of the triangle; if P is an arbitrary 

point of the triangle, then 

**(*> !/) = If) +«*&(!, If) + «3^ 3 (x, If) 



where 

_ area(PP 2 P 3 ) _ area (PP 3 P.) . _ area(PP,P 2 ) 
Vl area(flP 2 ft)' ' 2 aiealftftft)’ area(ftPrft) 



as shown in Figure 2.5.2. 
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P 

2 




Figure 2.5.2 - PiPjPi w split into three triangles 

Since <^3 are all positive when P is inside the triangle PiPtPs, all the 
coefficients in T, are positive: M« has the form 




(2.5.10) 



Now we assemble the element matrices to form the global mass and stiffness 
matrices. The membrane is replaced by an assembly of triangles A, with vertices 
Pi and edges P,P, as shown in Figure 2.5.3. The boundary condition u = 0 is 
imposed on the outer vertices labelled ‘O’. 




Figure 2.5.3 - .4n assembly of triangular elements 
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Figure 2.5.4 - The angles bettoeen out tow'd drawn normals to the faces are all 

obtuse 



We note that if i ^ j t then my >0, ki, < 0 iff P«. P> are the ends of an edge 
PiPj of the mesh. 

The finite element analysis of a 3D- acoustic cavity proceeds in a similar way. 
The elements are taken to be tetrahedra. and the pressure p(x . y, z) is taken as 

p|z, y, z) = a + bx + cy 4* dz (2.5. 12) 

in each tetrahedron. Now it is found Zhu (2000) [342]. (Jladwell and Zhu (2002) 
(131) that if the angles between the normals to the faces are all obtuse, as shown 
in Figure 2.5.1, then the element mass and stiffness matrices have the form 



This means that when the matrices are assembled they have the same kind of 
sign pattern as before: if t ^ j then niij > 0, kij < 0 iff P\Pj is an edge of the 



+ + + + 


" + 1 


+ + + + v 
+ + + + ’ K -~ 


1 1 
1 + 
+ 1 
1 1 


+ + + + 


. +J 
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Applying Lagrange's equation to the energies 

T = lu r Mu, V = iu r Ku 
2 2 

we find the equation governing the vibration as 

Mii + Ku = 0. (2.5.14) 

2.6 Natural frequencies and normal modes 

The matrix equation (2.4.13) represents a set of second order equations with 
constant coefficients. Following usual practice we seek the elution in the form 



q = 


9i ’ 
<h 




■ 

x 2 


sin(u* + tf), 


(2.6.1) 




. . 




In 

■ 







where the constants j? r . frequency u and phase angle $ are to be determined. 
When q has the form (2.6.1), then 

q = -^q = -u 2 xsin(o/t 4- £), (2.6.2) 

so that equation (2.4.13) demands that 

(K - AM)x = 0, \=J 2 . (2.6.3) 

This is the eigenvalue equation (1.4.1) and, since M is positive-definite and 
K is either positive semi-definite or positive-definite, the whole of the analysis 
developed in Section 1.4 can be used here. Thus the equation has n eigenvalues 
(A • )? satisfying 

0<A I SA 2 <...<A n! (2.6.4) 

and n corresponding eigenvectors (x,’)p satisfying 

(K - AiM)Xi = 0. (2.6.5) 

The frequencies cJi = (A«)t are called the natural frequencies of the system, and 
the eigenvectors are called the normal or principal modes. Note the distinction 
between Xi, a scalar, and x«, a vector. 

In order to become acquainted with the properties of natural frequencies and 
normal modes we shall consider the system specified by equation (2.2.7) and, to 
simplify the algebra, shall assume that 

( 2 . 6 . 6 ) 
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In this case the eigenvalue equation may be written 



where 



r 2 — A -1 0 0 




• • 

x\ 


O 

1 

* 

1 

1 




x 2 


0 ••• 2 — A -1 




Xn-1 


0 ••• -1 1 — A 

m 




*« 

. 



A = mu/ 3 /k. 



(2.6.7) 



( 2 . 6 . 8 ) 



To solve for the x T we use the idea suggested in Exercise 1.4.4, namely to write 
(2.6.7) as the recurrence relation 



-*r-i + (2 - A)* r - s r+l = 0. (r = 1,2 n). (2.6.9) 

The first of equations (2.G.7) may be written in this form if x 0 is taken to be 
zero; this may be interpreted as stating that the left hand mass (mo) is fixed- On 
the other hand, the last of equations (2.6.7) may be written in the form (2.6.9) 
if j„,i is taken to be equal to x u . Thus the end conditions for the recurrence 
(2.6.9) are 

x o =0 = x„ tl -x n . (2.6.10) 

The recurrence relation has the general solution 

x T = A cos r0 + B sin r0, (2.6.11) 



where, on substitution into (2.6.9) we find that 0 must satisfy’ 

cos(r - 1)6 + cos(r + 1)0 = 2ccs0cosr0 = (2 - A)cosr0, 
sin(r - 1)0 + sin(r + 1)0 = 2cos0sinr0 = (2 - A)sinr0, 



i.e., 

2 cos 0 = 2 — X. 

The end conditions will be satisfied if and only if 

>1 = 0 = sin(n + 1)0 - sinn0 = 2cos|(n + 1/2)0] sin 0/2, 
ao that the possible values of 0 are 



9 “ ft “T5TTr’ 

while the corresponding values of A are 
‘‘ = ’— 

Thus, in the «th mode, the displacement amplitude of the rth mass : 

_ . l(2, -l)nr, 



x r = sin r6, = 



(2n + 1) 



-)■ 



( 2 . 6 . 12 ) 



(2.6.13) 
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The modes for the case n = 4, which are shown in Figure 2.6.1. exhibit properties 
that are held by all eigenvectors of a tridiagonal matrix (such as that in (2.6.7) ) f 




For a proof of the convergence of this class of discrete models to the contin- 
uous beam, and for an estimate of the discretisation error on frequencies and 
mode shapes, see Davini (1996) |74|. 

(a) the ith mode crosses the axis (t — 1 ) times - the zeros at the ends are not 
counted; 

(b) the nodes (points where the mode crosses the axis) of the ith mode interlace 
those of the neighbouring ((t — l)th and (i -f 1 )t h) modes. 

If instead of being free at the right hand end, the system were pinned there, 
then the analysis would be unchanged except that the end conditions would be 

2-0 = 0 = x n . (2.6.14) 

In this case 0 would have to satisfy 



e = <S>, = ~, i = 1,2, ... ,n - 1 
n 



so that 




Chapter 2 



38 



and the corresponding eigenvalues, which we will label (A, )" , would be 




In the ith mode, the rth displacement amplitude is 



(2.6.15) 



yr = sm(r^) = sin( — ). (2.6.16) 

The two sets of eigenvalues (A,*)” and (A^)7~ 1 are related in a way which will be 
found to be general for problems of this type (see equation (2.9.10)), namely 

0 < A, < A? < - < An-1 < a£_, < A„. (2.6.17) 



Exercises 2.6 

1. Consider the beam system of Figure 2.3.1 in the case when (m,)", = 
m, <k,)i = k, f/,)o = t. Show that the recurrence relation linking the 
(tir)JJ* 2 may be written 

tv-2 - illr-1 + (6 - X)u r - 4u y + i +iv+ 2 = 0 

where A = muPP/k. Seek a solution of the recurrence relation of the form 

ti r = A cos rO 4- B sin r# -f C cosh rd + Dsinh 

and find 0, so that the end conditions ti_j = 0 = u 0 = u n _ 1 = 
are satisfied. Hence find the natural frequencies and normal modes of 
the system; i.e., a clamped-clamped beam. A physically more acceptable 
discrete approximation of a beam is considered in detail by (Jladwell ( 1962) 
[103] and Lindberg (1963) [215]. 



2.7 Principal coordinates and receptances 

Theorem 1.1.5 states that the vectors (xjj spa** the space of n- vectors, so that 
any arbitrary vector q(f) may be written 



q(0 = PlXl +P2X2 +-" + p n x n . (2.7.1) 

This may be condensed into the matrix equation 

q = Xp, (2.7.2) 

where X is the nxn matrix with the x, as its columns i.e., x, = , r»). 

The coordinates p\,p2,...,pr„ called the principal coordinates, will in general 
be functions of t: they indicate the extent to which the various eigenvectors x, 
participate in the vector q. The energies T, V take particularly simple forms 
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when q is expressed in terms of the principal coordinates. For equation (2.7.2) 
implies 





q = Xp. 


(2.7.3) 


that 


T = * (X$) T M(X*) = ip r (X T MX)f.. 


(2.7.4) 


But the element in row t, column j of the matrix X 7 MX is simply x ? Mxy and 
according to (1.4.12) this is zero if t ^ j . a, if t = j. Thus 

X 7 MX = dia«rfo 1 ,a a ,...,a B ), (2.7.5) 


s> that 


T = j{a,p? + a 2 ^ + -.a„p? 1 |. 


(2.7.6) 


Similarly 


V = ip r (X T KX)p 


(2.7.7) 


and 


X 7 KX = diap(A,a,, A 2 a 2 , . . . A^a,,) 


(2.7.8) 


that 


V = i{Aiajp? 4- A 2 a 2 ^ + - • - A „ar.p^}. 


(2.7.9) 



Equations (2.7.6), (2.7.9) show that the search for eigenvalues and eigenvectors 
for a symmetric matrix pair (M, K) is equivalent to the search for a coordinate 
transformation q — p which will simultaneously convert two quadratic forms 
q 7 Mq and q 7 Kq to sums of squares. 

We shall now use the principal coordinates to obtain the response of a sys- 
tem to sinusoidal forces. Equations (2.4.2) and (2.4.14) show that the equation 
governing the response to generalised forces (Q r )f is 

Mq + Kq = Q (2.7.10) 

where q = If the forces have frequency u and are all in phase, 

then Q and q may be written 

Q = 4-sin(ad + *), q = xsin(<«/f + ^). (2.7.11) 

In this case equations (2.6.1) - (2.6.2) yield 

(K — AM)x = *fr. (2.7.12) 

To solve this equation we express x in terms of the eigenvectors x,, so that 

x = *,x, +n 2 x 2 + ---!r„x r . = X»r, (2.7.13) 

where sri, it 2 , . . . , x, are the amplitudes of the principal coordinates pi,pz , . . . , p*. 
Substitute (2.7.13) into (2.7.12) and multiply the resulting equation by X T ; the 
result is 



X r (K - AM)X* = X r * = S. 



(2.7.14) 
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But now the matrix of coefficients of the set of N equations for the unknowns 
* 1 , tr 2 , . . . , ir B is diagonal, and the ith equation is simply 

(A, - A)o,ir, = E,, 

so that 

'< = M|br <2 - 7 - 

In order to interpret this result we consider the response to a single generalised 
force Q r . In this case 

Q= * = {0,0,...,*r,0,...,0), 3 = * r {x rl ,x*,...,Xn} 

*,Z,i 

x ‘~«.tA.-Ar 

and the sth displacement amplitude is: 

n 

X, = ^2 fl « Z « = Q r.4’r (2.7.16) 

»— 1 

where 

*•-£&*■ l2 - 7A7) 

The quantity a rj is the rettptonce Bishop and Johnson (19G0) (31] giving the 
amplitude of response of q M to a unit amplitude generalised force Q r . The fact 
that a rt is symmetric, i.e., 

<2.7.18) 

is a reflection of the reciprocal theorem which holds for forced harmonic excita- 
tion. 

Exercises 2.7 

1. Use the orthogonality of the (jqjy w.r.t the inertia matrix to show* that 

xfMq = pi<n. 



2.8 Rayleigh’s Principle 

(Consider a conservative system with generalised coordinates (^r)i vibrating with 
harmonic motion given by (2.G.1). Its kinetic and potential energies will be 

T = iq T Mq = J 2 coa 2 vtTo, 
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V = 



«jq r Kq =sin 2 u/tVo, 



where 



To = lx T Mx, Vo = ix T Kx. 



Since the system is conservative, 



( 2 . 8 . 1 ) 



T+V = const., 



so that 
and therefore 
This we may write as 



w 2 cos 2 wtTo + (1 — cas 2 to v t)Ki = const., 
u> 2 To = Vo. 



Vo _x‘Kx 



If the system is vibrating freely at frequency i o, then io must be one of the natural 
frequencies and x the corresponding eigenvector. If w = on, then A = A, . x = x, 

and T 

^ = ^ < 2S2 > 

which agrees with equation (1.4.5). 

Rayleigh's Principle states that the stationary t alues of the Rayleigh Quotient 



Xu = 



x Kx 
xHXx 



(2.8.3) 



uieuiftf as a function of the components (x r )J , occur when x is an eigenvector 
X|. The corresponding stationary tnlue of Xr is A,-. 

Proof. Rayleighs Principle has a long history - see for example Temple and 
Bickley (1933) (322) or Washizu (1982) (330). We shall state the proof in a 
number of ways because each is instructive. First consider Xr as a ratio of Vo 
and To and write down the partial derivative of this quotient w\r.t. x,. We 



dTo 

-j_- = m^x, + m r2 x 2 • 

SVo 






= krlXl + *»2Z2 



+ m rn z n , 
+ k n x,„ 



d Vo. 1 6Vo Vo dTo 1 \dVo 

K^7S>-3T8B:* 






7T\^v ~ Aw ^7/’ 



dTo\ 



so that, on inserting the expressions for 8Vo/9xr and dTo/dxr we obtain just 
the rth row of the matrix equation (2.6.3) with Xr in place of A. The complete 
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set of n equations which state that V 0 /7 o is stationary w.r.t. all the (x r )f is the 
matrix equation (2.6.3) which is satisfied when x is an eigenvector x* and A is 
the corresponding eigenvalue A*. 

Now express the energies in terms of principal coordinates. If 

pi = iTi sin(wt + 



then equations (2.7.6), (2.7.9) show that 

7o= ^{oin?+a2iri+-- + a n sS}. 



Vo = - {Aiamf + Aacfeir* + • - ■ + A, ,<»,,«* } . 

Since M is assumed to be positive definite, there is no loss in generality in taking 
each a* = 1 , then 

_ Arf+Agffg +—+*»«;■ 

R ~ • 

so that, in particular, 




(Aa - A,)sJ + ■ • • + (A„ - A,)t* 
*{ + *2 + ... + irj« 



(2.8.5) 



Since the A,- are labelled in increasing (or non-decreasing) order, the quantities 
A| — Ai , t = 2 , 3, • . . , n are non-negative, and so 



Ah > Ai. 

If Ai is strictly less than A 2 then equality occurs only when 1 T 2 = 0 = ... = 
x n , i e. f when the system is vibrating in its first principal mode. Equation 
(2.8.5) states the important property that whenever values are taken for (av)" , 
the values of the Rayleigh quotient will always be greater than Ai and (when 
Ai < A 2 ) will be equal to Ai only if the ratios x\ : Z 2 : - . . x» correspond to those 
of the first eigenvector xn, X2i,...x n i. Equation (2.8.5) show’s that Ai is the 
global minimum of X R . and it may be proved in an exactly similar way that 

Ah < A*, (2.8.6) 

so that A,, is the global maximum of \ R . 

If Xi is an intermediate eigenvalue, so that X\ < X* < A n , then 

. . _ - + £". 4+1 (a> - a*W 

* * 

In this case Ah will not be strictly less nor strictly greater than A » for variations 
of the xy. Ah has a saddle point in the ith mode (n j = 0, j £ t). However, for 
computational purposes it is important that the difference betw r een Ah and A, 
depends on the squares of the quantities irj. This means that if x is ‘nearly* 
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in the ith mode, so that the x } with j jt i are much smaller than :r,, i.e., 
*, « 1 = <Ke), then A„ - A, = 0(e a ). 

Since M is positive definite. x T Mx > 0. and the problem of finding the 
stationary* values of the Rayleigh quotient X R given by equation (2.8.3) is equiv- 
alent to finding the stationary* values of x r Kx subject to the restriction that 
x r Mx = 1. This in turn is equivalent to finding the stationary* values of 

F = x T Kx - Ax r Mx, (2.8.8) 

subject to x r Mx = 1. Here A acts as a Lagrange parameter. Note that 

dF n n 

— = 2J2 k r.z. - 2A TO r» I »> 

r *-l *-l 

&> that the set of equations dF/dx r = 0 yields equation (2.6.3), viz. 

(K — AM)x = 0. ■ 



2.9 Vibration under constraint 

The concept of a system vibrating under constraint is important in the solution 
of inverse problems. Suppose a sy r stem has generalised coordinates (g r )j\ but 
they are constrained to satisfy a relation 

/to l,qi,...,qn)=0. 

For small vibrations about q t = 0 = . . . = <&,, this relation may be replaced by 
q r d = <fi9i + d 7 q, + . . . + d n q„ = 0 , 



where 




In' 



Two of the most important constraints will correspond to a certain qr being 
zero, or two, q, and q,, being equal. Now suppose that the system is vibrating 
with frequency u, where tv 2 = A, and 



q = xsinut. 

Rayleigh's Principle states that the (natural frequencies) 2 will be the stationary 
values of F, given in equation (2.8.8) but now subject to the further constraint 

x r d = 0. (2.9.1) 



Thus we must find the stationary values of 



F = x r Kx - Ax r Mx - 2f'X r d. 



(2.9.2) 
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where v is another Lagrange parameter (the 2 is inserted purely for convenience). 
The equations dF/ttx, = 0 now yield 

Kx - AMx - nd = 0. (2.9.3) 

By comparing this with equation (2.7.12) we see that i-d is a generalised force; 
it is the force required to maintain the constraint (2.9.1). 

In order to analyse equation (2.9.3) we express x in terms of principal coor- 
dinates, using equation (2.7.13). Then 



KXn - AMX* - i'd = 0. (2.9.4) 

Multiply throughout by X r and use equations (2.7.5) and (2.7.8) which show 
that both X r MX and X r KX will be diagonal matrices; tbe rth row of the 
resulting equation is 



A^ir, - Aa r * r - I’b, = 0, r = l,2,...,n. 


(2.9.5) 


where 




b = X r d. 


(2.9.6) 


Equations (2.9.5) yield 




•*. 

' a,(Ar - A)’ 


(2.9.7) 


which, when substituted in the constraint (2.9.1); i.e., 




x r d = ir T X r d = ir T b = 0, 


(2.9.8) 



yields the frequency equation 

(2 “> 

The form of this equation has important consequences. Consider first the 
case in which none of the 6* is zero. The coefficients (h?/ot)J will all be pc&itive 
and the graph of B( A) against A will have the form shown in Figure 2.9.1. Since 
f?(A» + 0) is wry large negative, B(Ai+i — 0) is very large positive, and #(A) is 
steadily increasing bet wen A, and A*+i, £(A) will have just n — 1 zeros, (Aj)*" 1 , 
that interlace the A* in the sense that 



A* < a; < A,*!, (i = 1,2, . . . ,n - 1). (2.9.10) 

This inequality may be interpreted as follows: if a linear constraint is applied 
to a system, each natural frequency increases (or, more precisely, does not de- 
crease). but does not exceed the next natural frequency of the original system. 
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Figure 2.9.1 - The eigenvalues of a constrained system interlace the original 

eigenvalues 

If all the bi are non-zero then the inequalities in (2.9.10) are strictly obeyed. 
Now, however, suppose some of the 6, are zero; in particular consider the con- 
straint 

*i=0, (2.9.11) 

for which (6,) j = 0. In this case (**)£ are the principal coordinates of the 
system and the corresponding eigenvalues are 

AJ=A 1+It (i = l,2,...,n — 1). (2.9.12) 

If the constraint is 

*i = 0, 1 < j < n, 

then the principal coordinates are *i,*2< • • - n j+i> • - • >*n> ao that 

K=^> 1 = 1,2 = 1 —jtj + — 1 . 



If the constraint is (2.9.8) and a>me particular bj is zero, then equation (2.9.5) 
shows that 




1 = ) 



is a solution corresponding to A = A;. This means that a constraint (2.9.8) 
with bj = 0 does not affect the jth principal mode. Figure 2.9.2 shows the form 
of B(A) when 62 = 0- The graph may a) pass to the left of A2, in which case 
A'u < A 2' = A 2 ; or b) pass to the right, in whkh case A’^ = A 2 , A^ > X 2 . 

If two constraints are applied, then the constrained system will have n — 2 
eigenvalues (A")J~ 2 satisfying 



<*= 1 . 2 , 2 ), 
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where Aj are the eigenvalues of the system subject to one of the constraints. 
Thus 



K<K<K'<K+i<*i+ 2 , (t = l|2,... f n-2) 



or 

A< < Xi < A^ 2 , (i = 1, 2, . . . , n - 2). (2.9. 13) 




Figure 2.9.2 - The form of B( A) when 62 = 0; either a) Aj < A2, A2 = A2 or b) 

Ai = A2,A2 > A 2 . 



2.10 Iterative and independent definitions of 
eigenvalues 

In this section we take a closer look at the eigenvalues of (2.6.3) in relation to 
the Rayleigh Quotient 

A„ = (2.10.1) 

We assume that K is symmetric (it may or may not be positive semi definite) 
and that M is positive definite. The importance of the latter assumptions is 
that the denominator of (2.10.1) is never zero and always positive for all x ^ 0. 
First, we note that X R is a homogeneous function of x in the sense that 

Ar(oc) = A„(x), c* 0. 

This means that we can always scale any x sjo that the denominator of (2.10.1) 
is unity, i.e., 

x'Mx = 1. (2.10.2) 

The vectors x with this property constitute a closed and bounded subspace 
D, C V,,. Now consider the Rayleigh Quotient on Di; it is 

Aw = x T Kx. 



(2.10.3) 
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This is a continuous function of the variables x x , x 2 , . . . , x n on the closed bounded 
region D x so that, by W tier str ass f Theorem on continuous functions, it attains 
its minimum value on D Xy i.e., for some vector x € D x . (Recall the definition of 
a closed set S: if } is a convergent sequence in S then its limit lCt y f = y 
is also in 5.) There may be more than one such minimizing vector, but there is 
always at least one. which we denote by x x . The corresponding minimum value 
of X R we denote by X x . We have the result 

Ai = min x r Kx = x[Kxi. (2.10.4) 

Having found x x and Aj t we set up a new minimum problem: finding the mini- 
mum of x r Kx on the suhspace D 2 of D x that consists of vectors x orthogonal to 
x 1t i.e., x satisfying x T Mxj = 0. This subspace is again closed and bounded so 
that by Weieretiass’ Theorem there is a vector x 2 € D 2 which minimizes x r Kx 
on D 2 ; the minimum value is A 2 . We have 

A 2 = mm x r Kx = x 2 Kx 2 , (2.10.5) 

and x 2 Mx 2 = 1, xjMxi = 0. Since A 2 is the minimum of x r Kx on D 2 , a 
subspace of Di, A 2 cannot be less than Ai, i.e., A 2 > Ai. 

Proceeding in this way we find a set of vectors x, and numbers Ai,(i = 
1 , 2 , . . . , n) such that 

Xi = min x T Kx = x^Kx*, (2.10.6) 

j Wi. (2io7 > 

and A, 

This procedure is iterative: we cannot set up the minimizing problem that 
gives A 2 until we have found Xi , and generally we cannot set up the minimizing 
problem that gives A, until we have found x,,x 2 , There is another 
procedure in which we can find any A,,x, without first finding Xi.Xj,. . . ,x,_i; 
this is called the mdtfxndenl or minimal procedure. 

In the independent procedure we start as before: 

Ai = min x r Kx = xTKxi. 

xtD, 

Now” we return to the analysis of Section 2.9 relating to vibration under a con- 
straint. The inequality (2.9.10) shows that if none of the (&,)? is zero, then 
the first constrained eigenvalue, AJ, is strictly leas than A 2 . Equations (2.9.11), 

(2.9.12) show that if the constraint is xi = 0, then Af = A 2 . The quantity *1 
is the amplitude of the component of xi in x, and on premultiplying equation 

(2.7.13) by xf K we see that 

xfKx = TixfMxi = m. (2.10.8) 

Thus iti = 0 means that x is orthogonal to xi w.r.t. the matrix M; this is the 
constraint which yields the maximum value of Aj, namely A 2 . 
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Thus 

max min x r Kx = A, (2.10.9) 

» tSf 

where x J. d means x r Md = 0; the d which maximizes the minimum is xi. 

We may now extend this analyses to higher eigenvalues by using (2.9.13); 
thus 

max min x T Kx = As, 

di.dj »**>i 

and generally 

max min x T Kx = A*.j. (2.10.10) 

d..d* d, »di 

>xdi a. 

Again, the d‘s that maximize the minimum in the general case are d| = X| . d 2 = 
X9.---.dn_l =Xn_,. 

The minimax definition of eigenvalues seems to have been noted first by 
Fischer (1905) (88). The iterative and independent definitions of eigenvalues 
are discussed at length in Courant and Hilbert (1953) (61), and in the more 
specialised volume Gould (1966) (151). The motivation for Gould’s book was 
the search for lou<er bounds for eigenvalues; discretising methods like the finite 
element method almost always lead to upper bounds. 

Exercises 2.10 

1. Examine the arguments in Sections 2.9, 2.10 in the case when two eigen- 
values are equal, e.g., Ai = A 2 . 

2. Use the minimax procedure to show that if stiifhess is added to a system, 
i.e., the stiffiress matrix is changed from K to K', and x T K'x > x r Kx 
for all x € V,„ then none of the eigenvalues of the system decreases. Why 
can you prove this result only for Ai by using the iterative definition? 
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Jacobi Matrices 



Let »o one say that I have said nothing new; the arrangement of the subject is 

new. 

Pascal’s Penstes, 22 



3.1 Sturm sequences 

In this Chapter we will analyse the properties of the eigenvalues and eigenvec- 
tors of systems with the special tridiagonal mass and stiffness matrices met in 
Chapter 2. We will start by considering systems like that for the system in 
Figure 2.2.1, for which the mass matrix is diagonal and the stiffness matrix is 
tridiagonal, with negative codiagonal. At the end of the section we will show' 
that many of the results may be generalised to apply to systems like that in 
(2.1.10) in w'hich the mass matrix is tridiagonal with positive codiagonal. The 
most important property of the eigenvalues of such systems is that they are 
simple, i.e., distinct (Theorem 3.1.3). Thus 

Ai < Aa < ... < A*. 

If x r is the rth eigenvector, then as r increases, the eigenvectors oscillate more 
and more (Theorem 3.3.1) in such a way that the zeros of x r interlace those of the 
neighbouring and x r4 . 1 (3.3.4). We shall now establish these and other 
results. Throughout the next few Chapters, we redevelop analysis originally 
established by Cantmacher and Krein (19!>0) (98]. Their book was republished 

in 2002.. 

We start with a definition: 

Definition 3.1.1 A Jacobi matrix is a posifii'c semi-definite symmetric tridi- 
agonal matrix with (strictly) negative codiagonal. 

Note: Different authors define a Jacobi matrix in different ways; some choose 
the codiagonal to be strictly positive. 




so 
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Now we consider the equation 

(K - AM)x = 0 



(3-1.1) 



where K is a Jacobi matrix. First, we suppose that M is a (strictly) poritiw 
diagonal matrix, as in (2.2.7). and we reduce (3.1.1) to standard form. 

Take 

M = dxag{m it m 2 m„) 

and write M = D 2 , where 



D = diag(di.d2. ...,dn), d= mf. 



introduce the vector u related to x by 

u = Dx, x = D 1 u 



and premultiply (3.1.1) by D 1 to obtain 

D -1 (K - AD 2 )D -, u = 0, 



i.e., 

where 



(J - AI)u = 0, 



(3.1.2) 



J = D _, KD (3.1.3) 

The matrix J, like K, is a Jacobi matrix, and has the same eigenvalues as the 
system (3.1.1). We write 



ai — bi 0 ... 0 

— bi 02 —62 ... 0 



— 6»— 1 

-6»-i a» 



(3.1.4) 



The analysis now centres on the leading principal minora (see (1.4.6)) of the 
matrix J - AI. We define 

P 0 = l, P l (\) = a l -\ P 2 { A) = 

so that finally 

Pn( A) = det(J - AI). 

The minors satisfy the three-term recurrence relation 

P -+1 (A) = (or«i — A)P .(A) — b*P r-i(A), r = 1,2, .. . ,n — 1, (3.1.7) 

which enables us to calculate successively from Since 

the zeros of any P r ( A) are the eigenvalues of the truncated symmetric matrix 
obtained by retaining just the first r rows and columns of J, they are all real. 
We now prove: 



(3.1.6) 



V aj — A 
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Theorem 3.1.1 // > 0 (r = 1,2.... ,n-l), then the (P r (A»g form a Sturm 
sequence , the defining properties of which are 

1. P 0 (A) has everywhere the same sign (P 0 (A) = *)• 

2. When P r (A) vanishes, P r +i(A) and Pr- i(A) are non-zero and have opposite 
signs. 



Proof. In order to establish property 2 we note first that tw r o successive P r 
cannot be simultaneously zero - i.e., for the same A = A 0 . For if P*+i(A°) = 
0 = P,(A°) then equation (3.1.7) shows that P*_i(A°) = 0, so that finally w<e 
must have Pj and P 0 zero; but Po(A°) = 1, which yields a contradiction. 

The latter part of property 2 now follows directly from (3.1.7). ■ 

Before proceeding further we must define the sign change function s r (A). 
This is the integer-valued function equal to the cumulative number of sign- 
changes in the sequence P 0t Pi (A), P 2 (A), .... P r (A). Thus if 



J = 



2 - 10 ' 
-13-2, 
0-2 4 



then. 



Po = 1, Pt(A) = -A + 2, 

P 2 (A) = A 2 - 5A 4- 5, P 3 (A) = -A 3 4- 9A 2 - 21A + 12. 

For A = 0 the sequence of values is 1, 2, 5, 12. Since there is no change of sign 
in the sequence, each s r (0) = 0. For A = 3 the sequence is 1, -1, -1, 3, so that 
*i(3) = *a(3) = l, * 3 ( 3 ) = 2. 

Theorem 3.1.2 s r (A) changes only uhen A passes through a zero of the last 
polynomial, P r (A). 

Proof. Clearly. s r (A) can change only when A passes through a zero of one 
of the Pj(A), (s < r); it therefore suffices to prove that s r (A) does not change 
at all when A passes through a zero of an intermediate P,(A), (s < r). Suppose 
P,( A 0 ) = 0, where 1 < s < r, then P,_ i(A°) and P,*i(A°) will be both non-zero 
and have opposite signs. The signs of the triad P,_i(A°). P,(A°), P,+ i(A°) are 
therefore + 0 - or - 0+. Suppose the first to be the case, so that P # ( A) increases 
as A passes through A 0 , (the other possibility' may be handled similarly). Then 
for values of A sufficiently close to A 0 and less than A 0 the signs are 4* - 
while for values of A sufficiently clcee and greater than A 0 the signs are 4- 4- 
-. Thus, whether A is greater than or less than A 0 there is just one change 
of sign in the triad of values of P*-i(A), P.(A), P*+i(A). In other words the 
triad of polynomials P._i(A), Pi(A), Pj+i(A) will not contribute any change 
to s r (A) as A passes through A 0 . But no other members of the sequence will 
contribute any change to s r (A) as A passes through A 0 (unless A 0 is a zero of 
another Pi(A), \t — s\ > 2, in which case again there will be no change in s r (A)) 
ao that Sr (A) will not change at all. ■ 

Clearly, Sr(A) is not well defined when P r (A) = 0. 
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Theorem 3.1.3 The zeros of P y ( A), ore simple , L e. r distinct. In addition , if 
P r (A a ) / 0 and s r (\°) = A:, t/*en P r (A) has k zeros less than A 0 . 

Proof. Since P.(A) = (-)'A' H , all P.(A) will be positive for sufficiently 

large negative A. i.e., A < o t so that s r (o) = 0: a may be taken to be zero if 
J is positive definite. ()n the other hand, for sufficiently large positive A, i.e. t 
A > ft, the P*(A) will alternate in sign, so that &r(9) = r. Now since s r (A) 
can increase only when A passes through a zero of P r ( A), all the zeros of P r (A) 
must be distinct. For if A 0 were a zero of even multiplicity then s r (A) would 
not increase at all as A passed through A° f w'hile s r (A) would increase only by 
unity if A 0 were a zero of odd multiplicity. The second part of the theorem now' 
follows immediately. ■ 

Corollary 3.1.1 The eigenvalues of a Jacobi matrix are distinct. 

Corollary 3.1.2 The number of zeros of P r (A) satisfying o < A < (J is egual to 

* r (P) - *»•(<>)• 

Corollary 3.1.3 7/A° is a zero of P r (A) then , as A passes from A 0 — to A°+ 
the sign of P r -i(A)P r (A) changes from to and s r (A) increases by unity. 

Theorem 3.1.4 Between any two neighbouring zeros of P r (A) there lies one 
and only one zero of P r _i(A). and one and only one zero of P r +i(A). 

Proof. Let p^p^ ^ two neighbouring zeros. Suppose, for the sake of 
argument that P r (p x -) > 0. then P r (/ij+) < 0 and P r {p 2 -) < 0. By Corollary 
3.1.3, P r -iOi!+) > 0 and < 0, so that P r _!(A) change sign between 

pj+ and P 2 ~~i and therefore has at least one zero in (/i lf /i 2 ). 

Now’ property 2 of Sturm sequences shows that P r +\{Pi) and P r -i(p,)> (* = 
l t 2) have opposite signs. Thus Pr+i(p x +) < 0, P r4 .i(/i 2 -) > 0 so that P r + X ( A) 
has at least one zero in Now* suppose, if possible, that P r _i(A) (or 

P r+l (A)) had two (or more) zeros in (p x ,p 2 ) then P r (A) would have a zero in 
contrary to the hypothesis that p x ,p 2 are neighbouring zeros. ■ 

This theorem is usually stated in the form: the eigenvalues of successive 
principal minors interlace each other. 



3.2 Orthogonal polynomials 



There is an intimate connection between Jacobi matrices and orthogonal poly- 
nomials. In this section we outline some of the basic properties of orthogonal 
polynomials. 

Two polynomials p(x), q{x) are said to be orthogonal w.r.t. the weight 
function u*(x) > 0 over an interval (a, 6) if 




w(x)p(x)q(x)dx= 0 . 



(3.2.1) 
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A familiar example is provided by the Laguerre polynomials (L N (x))^ t i.e., 
Io(z) = 1, Li(x) = x — 1, L 2 (z) = x 2 - ix + 2, . . . 
which are orthogonal w.r.t. the might function over (0, oo), i.e. t 

J e-*Ln(x)Lr tl (z)dx = 0, m ^ n. 

One of the important properties of such polynomials is that they satisfy a three- 
term recurrence relation. The relation for the L r (x), for example, is 

Jt+i(*) = (*- 2 r - 1 )L r (x) - t^Lr.iix). 

In this section w'e shall be concerned, not with a continuous orthogonality 
relation of the form (3.2.1), but with a discrete orthogonality relation 

(p, q) = jr w,p{(.,m) = 0; («*0" > 0. (3.2.2) 

»-l 

where ({ g )f are n points, satisfying (i < < * * * < £*• 

To introduce the concept formally w r e let P„ denote the linear space of poly- 
nomials of order n, i.e., the set of all polynomials p(z) with degree k < n. with 
real coeflicients. On this space (,) acts as an inner product since it is positive 
definite , bilinear and symmetric , i.e., 

1. (p,p) = ||p|| 2 > 0 if p{x) ^ 0 

2. {ap,q) = a(p,q), {p + q,r) = (p,r) + (q,r) 

3. (p,9) = fa,p) 

In addition 

-i. (xp,q) = (p,zq) 

We now prove 

Theorem 3.2.1 There is a unique sequence of manic polynomials, Le. r (tfi(x))o 
such that Qi(z) has degree i and leading coefficient (of z*) unity, which are or- 
thogonal with respect to the inner product () r i.e., for which 

(<?.,&) = o. 3- 

Proof. The g g (x) may be constructed by appfying the familiar Gram- 
Schmidt orthogonalisation procedure to the linearly independent polynomials 
(xX’ 1 . Thus 

*— l 

«0 = 1 . *(*)=*' - (« = I|2,...,n — 1), 

>-o 
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where 

(«.9>) = (*\fc) = 0 

&> that 

«« = (**»«*)/! WP» J = 0, 1. ■ 

We note that the polynomial 

</..(*) = FI(z -{.) 

»- 1 

is the monic polynomial of degree n in the sequence. It is orthogonal to ($t>o"’ 1 ; 
in fact it is orthogonal to all functions, since = 0. t = 1, 2. . . . ,n. 

The (if am -Schmidt procedure does not provide a computationally convenient 
means for computing the instead we use Forsythe (1957) (90). 

Theorem 3.2.2 The monte polynomials (g*)© satisfy a three-term recurrence 
relation of the form 

q>(z) = [x - a*)$-i(x) - tf-iqi-2(x), i = 1,2,. ...n, (3.2.3) 

u nth the initial values 

q.,(x) = 0, qo(x) = l. (3.2.4) 

Proof. Qi( x )~ x ^i-i( x ) a polynomial of degree (a — 1). It may therefore be 
expressed in terms of (the linearly independent - see Ex. 3.2.1) 

Thus 

tt(*) - = < ‘o9o + c i9i + - • • + (3.2.5) 

Take the inner product of this equation with (j = 0, 1 ,. . . , i — 1) thus 

<9.- 1,) - (*-!.**) = 5Z c *(9*.9i) = Cjlfoll*. (3-2.6) 

k-a 

where the second term on the left has been rewritten by using property 4, above. 
But if j = 0,1,..., t - 1, then the first term on the left is zero, and if j = 
0, 1, . .. ,i — 3, then zqj has degree at nuret i — 2 and so is orthogonal to qt-i. 
Thus Cj = 0 if j = 0, 1,2, . . . ,i — 3 and there only two terms e,-i and c,_2 on 
the right of (3.2.5) i.e., 

9i(*) - **- 1(*) = ^-29.-2(^) + Ci_i9i-i(*). (3.2.7) 

Moreover equation (3.2.6) gives 

a* = -Ci-i = te-i,*»-i)/||flf-i|| a (3.2.8) 

Ci-2 = -( 9 . - 



i.*9<-2)/ll9.--2|| 2 . 
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But xq i ^ 2 is a monic polynomial of degree i — 1; it may therefore be expressed 
in the form 



• - 2 

*fc-2<*) = +'52 d ,<ij( x ) 

1-0 



so that 



(9.-1. -*9.-2) = Il9.-ill 2 

and thus negative and equal to — 3?_ lf where 



ft = IM/lta-ill- ■ <3.2.9) 



Equations (3.2.3), (3.2.4) with (3.2.8), (3.2.9) enable us to compute the poly- 
nomials {^t)7“ l step by step. Thus with q - 1 = 0, <jo = l we first compute ai 
from (3.2.8); this substituted into (3.2.3) gives <71. Now we compute 0-2, ft and 
find </2, etc. 

In inverse problems we will need to express the weights t l\ in terms of the 
polynomials q n _ x and q n . For this we note that if f(x) is any polynomial in 
P n _i, i.e., of degree n — 2 or less, then 

<9n-l./) 3 jrw.qn-liQfiti) = 0. 
i-l 

But if such a combinat ion 



"*• = «®i9a-i(0 

• -1 



is zero for any f(x) in P*_i then 

n 

= (fc = 0 ,l,...,n — 2 ), 



t-1 



x k is in P n— I. i.e., 



1 1 1-1 

Bm=| il 

q-2 ^n-2 { n-2 ... 

U is shown in Ex. 3.2.2 that this equation has the solution 

-= 7 /n'(^,). 





mi 




m 2 




ms 




m n 



= 0 . 



(3.2.10) 



(3.2.11) 
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Apart from the arbitrariness of the factor 7, this is the unique solution. The 
prime means that the term j = i is omitted. Now since 

9,(0 =fl «-<,). (3-2.12) 

7-1 

we have 

<.cu -n 

j-i 

where the prime on the left denotes differentiation! 

Returning to equation (3.2.11) we can deduce that, for some 7, 

«*. = f.9,-i(0) = W,(0)- * = *- 2 n - 

Since the {$}£ satisfy the three-term recurrence relation (3.2.3) it follows, 
by the arguments used in Section 3.1, that the zeros of ^ n (x) and q,, i(x) must 
interlace and therefore (Ex. 3.2.3) 9„_i(£,)9i(ij) > 0. This means that the 
weights 

** = -y/fo«-i (&)£(&» <3.2.13) 

are positive. 

This equation is important: it means that if the monte polynomials 

?„((), q n -i(£) ar * given, and if their zeros interlace , then they may he viewed as 
the nth and (n — l)th members , respectively , of a sequence of monic polynomials 
orthogonal w.r.t. the weights us given by (3.2.13), and the points {(<}*. 



1. Show that if the polynomials {£{}§, k < n , are orthogonal w.r.t. the inner- 
product (3.2.2). then they are linearly independent. Hence deduce that 
any polynomial p(z) of degree A: — 1 may be expressed uniquely in the form 

it- 1 

p|i) = 53c i9i (i), 

7-0 

and that qk(x) is orthogonal to each polynomial of degree k — 1 . 

2. Show that if the Vandermonde determinant A is defined by 





1 


i 


i 




(i 


& ■ 


• £..-i 


A = 


d 


d • 


■ d-i 






a~ 2 • 


• «:? 
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A="ri n«,-«.>=r/n, 

>-2 t-l )-l 

r =ii 

7-2 fc-1 



Hence deduce (3.2.11). 



3. The zeros {&}f of g«(x), and {&}" -1 of <fr-i(r) must satisfy (i<(I< 
& < - < C-i < <«• Show that (-)-VJ^) > o, iK() > o 

and hence 9 ' ri K,) 9n _ 1 ({.)>0. 



3.3 Eigenvectors of Jacobi matrices 

In this section we establish some properties of the eigenvectors of Jacobi matrices, 
in preparation for the solution of ‘inverse mode problems’. We return to the 
analysis of Section 3.1 and prove 

Theorem 3.3.1 The sequence (u rj -)p_ 1 for the j th eigenvector has exactly j— 1 
srgn reversals. 

Proof. The u,j are determined bom equation (3.1.2) for A = Aj; this may 
be written 

-hr-r^r-ij + (a. - “ Mr+ij = <r= 1,2,..., n) (3.3.1) 

where uo,>,Un«i,; are interpreted as zero, i.e., 

tioj = 0 = i*n*ij. (3.3.2) 

Choose an arbitrary 6 n > 0 and put 

t'r = **iji W = hiuaj,..., ff.«I = 6|6j • hatin.l J 

and multiply equation (3.3.1) by b t bj ■ • ■ b,_ l to obtain 

+ (a r -A^K-t* P4I =0, (r = 1,2, ...,n). (3.3.3) 

()n comparing this equation with (3.1.7) we see that it has the solution 

"0=0, IH = 1, W = Pr-i(hi), ( r = 1,2, ...,n+ 1) 

which because of P n (Aj) = 0, satisfies the end-condition Vn+\ = 0. 

Thus, 

Urj = (6lk ••6 r _l)- , Pr-,(A i ), (3.3.4) 

and since A* lies between the (j-l)th and jth zeros of P n -x (\), = j-1. 

■ 

Before establishing further properties of the eigenvectors we introduce the 
concept of a a-line. 
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Definition 3.3.1 Let u = {u 1 ,u a ,...,ti a + l } he a vector. We shall define the 
u-line as the broken line in the plane Joining the points with coordinates 

x r = r, y, = u„ (r = 1,2 ,... ,n + 1). 

Thus, between (x„y,) and (* r+1 ,!' r +i), v(*) * defined by 

y{x) = (r + 1 - i)u r + (* - r)Ur+ t , (r = 1,2, . .. ,n), 
as shown in Figure 3.3.1. 

Now return to Theorem 3.3.1. For arbitrary (real) A the sequence given by 

no = 0, u,(A) = u,(b,b 3 ■ ■ • 6r_,)- , / > r _,(A), (r = 1,2, . .. ,n + 1), 

satisfies the recurrence (3.3.1) for r = 1,2 (It will satisfy the last 
equation with tq*+ x = 0 iff P n ( A) = 0.) For arbitrary' A, the vector u(A) = 
(tii(A), . . . ,tq 4 + i(A)} defines a u(A)— line. We now' investigate the nodes of this 
line, i.e., the points x at which y(x) = 0. First we note that if ti r (A) = 0, i.e. t 
/’r-i(A) = 0, then P r ( A) and P r . 3 ( A), i.e., u rtl and «•_, will have opposite 
signs, so that the ti(A)— line will cross the x-axis at x = r. Secondly, if tq. 
and tq.+ x haw opposite signs, then y(x) has a node between r and r 4- 1. This 
implies that the u(Aj)— line has exactly j nodes, excluding the left hand end 
where no = 0, but ineluding the right hand end. Moreover, if A/ < A < Aj+i, 
then the u(A)— line will have exactly j nodes, again excluding the left hand end 
where no = 0. Table 3.3.1 shows the signs of Ur for the whole range of A- values, 
for the case n = 3. The last line in the table shows the number of nodes in the 
ti(A). Figure 3.3.1 shows the form of the ti(A) for the starred values of A. We 
now establish an identity which will enable us to prove further results concerning 
the eigenvectors. 
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Table 3.3.1 -The signs of u, /or different lvalues of A 
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Ai 
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A 


Aa 


A” 
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+ 


+ 
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+ 


0 
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0 
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0 


+ 


**4 


+ 


0 


- 


- 


0 


+ 


0 




0 


1 


1 


1 


2 


2 


3 



Consider the solutions u,v of the equations (3.3.1) corresponding to A, p 
respectively. Suppose that uo = 0 = no and that some positive value has been 
aligned to b,,. Then 

-6 r _iu r _ I +a r u r -6 r u r + I = Au,, (r = l,2,...,n) 

-6 r -i*V-i +o,u r -6 r ti r+ i = itv r , (r = l t 2,...,n). 

Eliminating a, from these equations, we find 

Mr - Mlf r -l = </l - A)(«r - «r-l) 

where 



»r = - tlrtV+1, Sr = 

i-1 



(3.3.5) 

(3.3.6) 



so that on summing over r = p,p + 1 ,. . . .q (1 < p < q < n), we obtain 

M* - = (P~ A H S « - «)- (3-3-7) 

In particular, if p = 1, so that Uq = 0 = Vq = t 0 = *o> 

V, = 0-A)s,. (3.3.8) 

We now prow 



Theorem 3.3.2 If A < /i, then between any two nodes of the u(A) - tine there 
is at least one node of the u(fi) - line. 

Proof. Let a t (3(o < ft) be two neighbouring nodes of the u-line and suppose 
that 

p-l<a<p, q<3<q+ 1, (p<<?), 

so that 

y{a, A) = (p- a)uy>-i + (a - p + l)ti* = 0, (3.3.9) 

y(A AJsh+l-flK + W- q)«i+i = 0 (3.3.10) 

and y(x,A) ^ 0 for a < x < 0. For the sake of definiteness suppose that 
y(z, A) > 0 for a < x < ft y then u Vl tip+i, . . . . ti v are all positive. We now need 
to prove that y(x,/i) has a zero between a and ft- Suppose y(z,/i) has no such 
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zero, that is, it has the same sign for a < x < 0 . Without Id's of generality we 
can assume that 

y(z, p) > 0 for o < x < 0 . 

that is y{ct,p) > O,p(0,p) >0 and tyn,. ..,v v are all positive. Thus. 

(p- o)*v-i + (a-p+l)up>0, (3.3.11) 

(<? + 1 - 0 )v q + [0 - q) tv+i > 0 , (3.3.12) 

and on eliminating a bet wen (3.3.9), (3.3.11). and & between (3.3.10), (3.3.12) 
we deduce that t v -\ > 0. tq < 0. On the other hand, s^-Sp-i = W > 0, 
so that the LHS of (3.3.7) is non-positive, while the RHS is positive, providing 
a contradiction. If we had assumed y(x./i) < 0 for a < x < the we wx>uld 
haw found the LHS of (3.3.5) non-negative and the RHS negative. ■ 

Theorem 3.3.3 As A increases continuously . then the nodes of the u(A) - line 
shift continuously to the left. 

Proof. Let o l (A),a 2 (A),. .. be the nodes of the u(A) - line, and suppose 
0 < (k\ (ji) t a 2 (/i), . . . are the nodes of the u(p) - line. We need to prove that 

a r (/i) < o r (A) 

for all these values of r corresponding to the u(A) - line. Since, by Theorem 
3.3.2, there is a least one of the o r (/i) between any two of the a r (A), it is sufficient 
to prove that 

ai(/i) < ai(A) = x . 

Suppose if possible that &i(p) > x and that 

q<x<q + 1 (1 < 9 < n) 

then all . . , i* v and . , u, will be positive while 

(9+l-*)«, + (*-9)u, +l = 0 
(9+1 -*)», + (*-9)i», + , > 0 

which imply t, < 0. On the other hand a, > 0, which, when used with (3.3.8), 
provides a contradiction. ■ 

Table 3.3.1 shows how the first node of ti(A) appears at the right hand end 
(n+ 1) when A = Ai and gradually shifts to the left, how the second zero appears 
when A = A2, etc. 

Theorem 3.3.4 The nodes of two successive eigenvectors interlace. 

Proof. Let the eigenvectors correspond to A, and A,,i. The nodes of the 
«*(Aj) and u(A,,i) - lines are (<tr(A,))'_, and (o r (A„i))Jij respectively; and 
= Qj*.i(Aj*.i) = n + 1. Theorem 3.3.3 shows that oi(A^*i) < ai(Aj), 
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while Theorem 3.3.2 applied to the two zeros a^.x(A^) and a^(A^) = n+ 1 shows 
that These two inequalities imply that the only possible 

ordering of the nodes is 

0 < ^(A^) <a 1 (A i ) <a a (A i+1 ) <••• < Qj-riAj) (3.3.13) 

< < MV) = a*+i(V+0 = n + 1 . ■ 

The derivation of certain other important properties of the eigenmodes will 
be deferred until Section 5.7, where properties of an oscillatory matrix will be 
used. See Gladwell (1991a) [119] for some related results. 

Exercises 3.3 

1. Show that the first and last components of any eigenvector of a Jacobi 
matrix must be non-zero. 

2. Show that if the matrix J of (3.1.4), with negative off-diagonal elements 
has an eigenpair A,,u,\ then the corresponding matrix J # with positive 
off-diagonal elements, has eigenpair A*, Zu* where Z is given by Z = dxog 
(1, —1, 1, • • - (— )"“ l ). This means that the eigenvector corresponding to 
the smallest eigenvalue, Ai , has n— 1 sign changes, w'hile that corresponding 
to A,» has none. Show* that if the eigenvalues of J* are numbered in reverse, 
i.e., Ai > A 2 > • •• > A,, > 0, then Theorem 3.3.1 remains valid. 

3.4 Generalised eigenvalue problems 

In Section 2.4 we showed that the eigenvalue problem for a finite element model 
of a vibrating rod could be reduced to a generalised eigenvalue problem 

(K — AM)u = 0 (3.4.1) 

where K,M were both symmetric tridiagonal matrices. K having negative co- 
diagonal and M having positive codiagonal. If M is positive definite, and K is 
positive semi-definite (i.e., K is a Jacobi matrix), then the analysis of Chapter 
1 shows that the eigenvalues are non-negative. Under these conditions we may 
prove that the solutions of (3.4.1) share the properties of the eigenvalue prob- 
lem in normal form i.e., equation (3.1.7). In particular, we can show that the 
eigenvalues of (3.4.1) are distinct and that the sequence (ti rj )p_ 1 for the jth 
eigenvector has exactly j — 1 sign reversals. To obtain these results we need to 
return to the analysis in Section 3.1 onwards and see what changes have to be 
made. 

We start with the principal minors of the matrix K - AM, using the notation 
of (2.4.10): 

ci — Aoi — di — A6i 

— dl — Atl C2 — Ad2 ’ ’ " 

(3,1.2) 



ft(A) = l, ft(A)=ci-Aai, ft(A) = 
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so that finally 



Pn( A)=det(K-AM). 



The minora satisfy the three-term recurrence relation 



Pr. 1 <A) = («V + 1 - Ao r «l)P r (AJ - (dr + A6,) J Pr-l(A). 



(3,1.3) 



The argument used in Section 3.1. 3.3 was based on the fact that the sequence of 
principal minors defined by (3.1.5) t (3.1.7) was a Sturm sequence. The sequence 
defined by (3.4.2), (3.1.3) however, is not a Sturm sequence. For if P r (A) = 0, 
then P r +i(A) = — (A - f A&r) 2 P r -i(A), and if it happens that d r 4- A A = 0, 
then P r +i(A) would be zero, and not. as required by condition 2 of Theorem 
3.1.1, of opposite sign to P r _i(A). Now we make the crucial observation, that 
if we restrict attention to A > 0, then the P r (A) do form a Sturm sequence 
because 6 r . A being positive, eliminates the possibility that d r + A b r = 0. If 
we assume that M is positive definite and K is positive semi-definite, then all 
the eigenvalues Ay will be non-negative and we may proceed as before. Thus 
Theorem 3.1.1 holds provided that A > 0. and Theorem 3.1.2 holds. The proof 
of Theorem 3.1.3 must be slightly changed. In the expansion of P,(A) in powers 
of A we have 

/’.(A) = Q,.o + <*.,iA + . . . + a,..A'. (3,1.4) 

The first term, a,.o, is the sth principal minor of K and, since K is positive semi- 
definite, q,,o > 0 for s = 1, 2,. . . , n — 1 and a B| o > 0; since P.(0) = o,.o we haw 
«,(0) = 0. The last term in (3,1.4), is a,., = (-)' * (the sth principal minor of 
M), so that for sufficiently large A, i.e., A > 0, Sr(0) = r. The remainder of the 
proof of Theorem 3.1.1, the corollaries 1-3 and Theorem 3.1.4 follow as before. 

We need to make small changes in the proof of Theorem 3.3.1. The iq.^ are 
determined from the equations 

+ Ajh,.,)**,.-!^ + (<V - A jOr)Ur.i - K + A^K^j = 0 (3.4.5) 

for r = 1,2, . . . ,n, where = 0 = u„ +M . 

Put d T + Xjb, = e„ choose an arbitrary e n > 0, and put 



«<2 = « n ., =e 1 e 2 ...e„ij n>lj 

and multiply equation (3,15) by eie 2 . . . e r _, to obtain 

+ (<V - AjO r )t' r - u r+1 =0 r = 1,2 n. 

On comparing this with (3.4.3) we see that it has the solution 

t\> = 0,n = l,tv = Pr-i(A>),(r = 1,2,... ,n + 1). 

Again, we conclude that *n-i(Aj) = j — 1. 

We may make similar changes to the proofs of Theorems 3.3.2-3.3.4. 



Exercises 3.4 



1. Make appropriate changes in the proofs of Theorems 3.3.2-3.3.1. 




Chapter 4 

Inverse Problems for Jacobi 
Systems 



People are generally better persuaded by the reasons which they themselves 
have discovered than by those which have come into the minds of otheis. 

Pascal's Penates, 10 



4.1 Introduction 

Research on these inverse problems began in the former Soviet Union, with the 
work of M.G. Krein. It appears that his primary interest was in the qualitative 
properties of the solutions of, and the inverse problems for, the Sturm-Liouville 
equation (see Chapter 10), and the discrete problems were studied because such 
problems were met in any approximate analysis of Sturm-Liouville problems. 
Krein’s early papers Krein (1933) (198), Krein (1931) [199] concern the theory of 
Sturm sequences, while the Supplement to Gantmacher and Krein (1950) (98), 
Gantmacher and Krein (2002) and Krein (1952) (202) make use of the theory of 
continued fractions developed by Stieltjes (1918) (310). Krein sees his results as 
giving mechanical interpretations of Stieltjes’ analysis. 

(Consider the simple system shown in Figure 4.1.1a. 




Figure 4-11 * Th* system is a) free and b) fired at the right hand end 

If mi, m 2 , k\,k2 are given, then the analysis of Chapter 2 shows how we can 
find the two natural frequencies, of the system: Ai = A 2 = are the 
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eigenvalues of the equation 

[ h+ V“ ][;;]=“ w 

The eigenvalues are the roots of the determinant: 

A(A) = mimaA 2 - {fami 4- (ki 4- Ar 2 )m 2 }A 4- k\k2 = 0. (4.1.2) 

Now consider the inverse problem. First it is clear that if one 9et of values 
wi,m 2 , Ai, has been found that yield specified eigenvalues Ai,A 2 , and if a > 
0, then ami,am2,aAi,a& 2 will be another set yielding the same eigenvalues: 
there are not four quantities to be found, only three ratios mi : m 2 : fa : fa. 
Knowing these ratios, we would need one more quantity, for instance the total 
mass m = mi 4- m 2? or the total stiffness k given by 1/A: = 1/fa -f 1/fa, to find 
the absolute values of the four quantities mi, m 2 , fa, fa. 

But even knowing two eigenvalues Ai, A 2 , wo cannot find the three ratios; wo 
need one more piece of information. One possible piece is the single eigenvalue 
A* = of the system obtained by fixing m 2 , as showTi in Figure 4.1.1b. This 

A»=o/ = A:i4 '*r (4.1.3) 



The sum and product of the roots Ai, A2 of equation (4.1.2) are 



Ai -f A 2 = 


A-ami + (&i + *T2)m2 *2 ( (fci + *2) 


(4.1.4) 




m 1 m 2 m^ 




, , *>** 
A1A2 = . 

mim 2 

Subtracting (4.1.3) from (4.1.4) we obtain 


(4.1.5) 




— = Ai + Aj - A*, 
m 2 


(4.1.6) 


and then (4.1.5) gives 
and finally (4.1.3) gives 


fci _ A1A2 
mi Aj 4 - A 2 — A 


(4.1.7) 



* 2 _ _ A* _ _ a* _ *1*2 = (A’ ~ A,)(Aa — A*) 

mi mi Ai + A2 - A' Ai + A2 - A' 



(4.1.8) 



The general theory of vibration under constraint (Section 2.9) states that Ai < 
A* < A 2 , so that all the quantities on the right hand sides of (4.1.6>-(4.1.8) are 
positive: the solution is realistic. The theory presented in this Chapter provides 
various generalisations of this analysis to a lumped-mass system made up of n 
masses. The Chapter falls into three parts: a discussion of inverse problems for 
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a Jacobi matrix; masa-spring realisations of these problems; generalisations and 
variants of these problems. 



1. Show that if ui,u 2 are the eigenvectors of (4.1.1), normalised so that 
ufMuj = Sij 9 then the equation giving the eigenvalue A* is 



** 2.1 



= 0 



Ai - A A 2 - A 
so that knowing A* is equivalent to knowing tr 2 . 2 : ti 2 ,i. 

2. Show that for the system of Figure 4.1.1, the system of given stiffness k , 
(1/A* = l/ki + 1 /fc 2 ), and least mass m = m r f m 2 , is found for 

A # = Aj -4- A 2 — v/^i^2* 

3. Show that for a taut string with tension T and unit length with just one 
concentrated mass m located at a distance t x from the left hand end. t 2 
from the right, the frequency u is given by 

k\ 4- k 2 — Am = 0 




6 + 6 = 1 . 



Hence find the system of least mass ha\ing a given frequency u = >/X7 
This suggests the problem of finding a string of least mass having concen- 
trated masses (mjy separated by distances ti f i 2f . . . ,£*+ 1 , where Ya-i ^ 
= 1. Barcilon and Tfrrchetti (1980) [23) considered this problem in a wider 
context, but did not find a closed form solution for the discrete problem. 



4.2 An inverse problem for a Jacobi matrix 

It was shown in Section 3.1 that the (natural frequencies) 2 of a lumped mass 
system may be obtained as the eigenvalues of a Jacobi matrix 



01 -61 

—61 02 -t>2 

0 0 On-1 -bn-l 

0 0 —bn-l a,. 



(4.2.1) 



If the system is connected, i.e., the stillnesses between mass 
tive, then the codiagonal elements -6,- are strictly negative. 
The basic theorem is 



are strictly posi- 
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Theorem 4.2.1 There is a unique Jacobi matrix J having specified eigenvalues 
(*)?, where 

0< Aj <A 2 <*'*<A n (4.2.2) 

and with normalised eigenvectors having non-zero specified tt dues or 

(t*m)i of their first or last components respectively; recall that u, = . . . , 

“»<>• 

(We recall Ex. 3.3.1, that the first and last components of an eigenvector of 
a Jacobi matrix are both non-zero.) 

Proof. The theorem is at once an existence (there is ...) and a uniqueness 
(... a unique) theorem. We shall prow existence by actually constructing a 
matrix, and will do so by using the so-called Lanraos algorithm; the algorithm 
demonstrates that J is unique. This algorithm has the advantage that numeri- 
cally it is well conditioned. An independent proof that the matrix is unique is 
left to Ex. 4.2.2. The proof will be presented for the case in which (u^)" are 
specified. 

The eigenvectors u« satisfy 

Ju. = AiU. (4.2.3) 

Use the column vectors (Ui) J to construct a square matrix U : U = (ui , 112 . . . . , u,,] 
The orthonormality conditions u^u> = Sij yield 

U T U = I. 



This means that U T is the inwrse of U : U is an orthogonal matrix. But if 
U T U = I. then Theorem 1.3.6 states that UU r = I also. Now put U r = X, 
then Ulr = X T X = I. But this means that the columns of X. like the columns 
of U, are orthonormal. Call the columns ( Xi)\ , so that X = (xi,X2, . .. ,x„), 
then 

xf X, = Sij. 

The reason why we have introduced the vectors x, is that 

xi = = {1111,1m,. ( 4 . 2 . 4 ) 



is given as part of the data. 

Now we proceed to rewrite the eigenvalue equations (4.2.3) as equations for 
the x,. The set of equations (4.2.3) for t = 1,2, .. . ,n, may be written 

JU = UA. (4.2.5) 



Thus, on transposing we find 

XJ = AX. 

Written in full, this equation is 

a, -6, 

-6, a 3 —6a 

-*n-i «*,. 



(4.2.6) 



(x,,x 2 ,... l x T1 ] 



= A(x ll xa,...,xJ. (4.2.7) 
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Take this equation column by column. The first column is 

«ixi — 61X2 = Axi. (4.2.8) 

Premultiply this by xf, using xfxi = l t xfx2 = 0; 

a,xfx, =a, =xfAx,. 

Now rewrite equation (4.2.8) as 

61x2 = mxi - Axi = Z2. 

The vector Z2 is known, because ai, xi , A are all known. The vector xj is to be 
a unit vector, so that 

M|x a || = 6, = ||z,|l- 

and X2 = za/fci. Hating found ai.&i,X2 we proceed to the next column of 
(4.2.7): 

-6lXl + 02 X 2 - 62X3 = AX 2 

Again, premultiplying by xj we find a 2 = x|" Ax 2 , and then 
62X3 = 02 X 2 - 61X1 - AX2 = Z3 

*> that 

&2IIX3II = &2 = ll>sll. *3 = Z3/63 

and so on. This procedure is called the La n c2QS algorithm; see Lanczos (1960) 
(203], Golub (1973) (132], Golub and Van Loan (1983) (135] and Kautsky and 
Golub (1983) (192). It produces a matrix J and at the same time constructs 
the columns (Xi)J which yield X = U T . 

Actually, what we have described is an inverse version of the original Lanczos 
algorithm. This original algorithm solved the following problem: Given a 

symmetric matrix A and a vector xi such that xfxi = 1, compute a symmetric 
tridiagonal matrix J and an orthogonal matrix X = (xi,X2 , . . . ,x«) such that 
A = XJX r . In our use of the algorithm, we start with A = A. 

We haw defined a Jacobi matrix as a positive semi-definite symmetric tridi- 
agonal matrix with strictly negative codiagonal. If the spectrum (Ai)" satisfies 
the inequalities (4.2.2), so that Ai > 0, then the J constructed by the Lanczos 
algorithm from A = A will be a Jacobi matrix. ■ 

Exercises 4.2 

1. Show that the vectors x, constructed in the Lanczos algorithm satisfy 

xTx i =6, i i,j = 1,2, ...,n 

even though this orthogonality is apparently established only for |t— j\ < 1. 
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2. Show that there cannot be two distinct Jacobi matrices J and J' with 
<r(J) = <y(J r ) and with the same values of the first components (u lt )f of 
their normalised eigenvectors. 

3. Rewrite the procedure described in equation (4.2.5) on, to solve the original 
Lanczos problem. 



4.3 Variants of the inverse problem for a Jacobi 
matrix 

First, we introduce some notation. Suppose A € M,*. The set of eigenvalues 
of A, the spectrum of A is denoted by o (A). If A is symmetric, i.e., A € S n , 
then tf(A) is a sequence of real numbere (AJf. where Xi < A 2 < A 3 • • • < A„. If 
M,K € then the set of eigenvalue of equation (3.1.1) is denoted by <r(M,K); 

again it is a sequence of real numbers (A*)? satisfying Ai < A 2 < • • • < A,». 

See Kautsky and Golub (1983) [192], deBoor and Saff (1986) [76] for a dis- 
cussion that places the Jacobi matrix problem in a wider context. FYiedland 
and Melkman (1979) [94] discuss the inverse eigenvalue problem in the context 
of non-negative matrices. 

If A € S ny the matrix obtained by deleting the tth row’ and column of A is 
called a truncated matrix. It will sometimes be denoted by A,; its eigenvalues 
will be denoted o(Ai). 

Now* suppose that A € S n is a Jacobi matrix J, then its eigenvalues will be 
distinct, and the eigenvalues <r(Ji) = (pi)\ 1 will strictly interlace (A,)", i.e., 

0<A I </i 1 <A 2 <-.-</.,,_ I <A n . (4.3.1) 

The problem of reconstructing J bom a( J) and <r(J , ) seems to have been studied 
first by Hochstadt (1967) (173). He proved that there is at most one matrix J 
with the required property. Hochstadt (1973) [176] attempted to construct this 
unique Jacobi matrix, but he did not show that his method would always lead to 
real values of the codiagonal elements b,. Hald (1976) [160] presented another 
construction and showed that, in theory, it would always work provided that the 
eigenvalues satisfied the interlacing condition (4.3.1). In practice, however, the 
construction was found to break down due to loss of significant figures. Hald also 
showed that Hochstadt’s construction will always lead to real b, provided that 
(4.3.1) holds. Gray and Wilson (1976) [154] presented an alternative, inductive 
construction of J. An independent uniqueness proof was given by Hald (1976) 
(160). 

In this section we shall present two methods for constructing J. The first 
relies on the theory of orthogonal polynomials described in Section 3.2. The 
second, which will later be generalised to inverse problems for band matrices, 
relies on the Lanczos algorithm described in Section 4.2. 

Note that we have chosen to define a Jacobi matrix so that it is positive 
semi-definite. Many of the results require only the interlacing of the A T s and 
the p’s, without any restriction on positivity. 
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The first method is best described by supposing that o( J) = (A t *)y and 
<y(J n ) = (/i,)?" 1 are known, and satisfy (4.3.1). Remember that J n is obtained 
by deleting the nth row and column of J. Now we form monic polynomials 
p t (A) t rather than the polynomials P*(A) in equation (3.1.5). We form two 
polynomials 



Pn(A) = n(A-<U Pn-.(A)=n< A -"<>- 



(4.3.2) 



t-1 



t-1 



The polynomials are the nth and (n — l)th monic polynomials of the sequence 
of monic polynomials with weights given by equation (3.2.13), i.e. f 



•".^/{P-iW^A,)} (4.3.3) 



and points (A,)". In addition, they are the nth and (n — l)th principal minors 
of the matrix (AI — J). The polynomials pr( A) therefore satisfy 

Pr(A) = (A - Ur)jV-l(A) - ^_ 1 p r - 2 <A). (4.3.4) 



Hald’s method of reconstructing J is as follows: he starts from p n (A).p PI _ 1 (A) 
and constructs p M - 2 (^)t ^d in the process finds a n and 6 n _j, by synthetic di- 
vision. Then from p K _i(A),p n _ 2 (A) he constructs p n _3(A) and finds a*_i and 
&*~2t “d so on. The process is inherently unstable because the polynomials 
Pn-2»P«-3i • •• >Pi found by successively cancelling the leading terms in the 

preceding pair of polynomials; the process becomes unstable because of cancel- 
lation of leading digits. 

de Boor and Golub (1978) [75] proceed quite differently. Having found the 
weights by using (4.3.3), they construct the polynomials in the natural order 
by using the analysis of Section 3.2, i.e., 

P-i(A) = 0, po(A) = 1 (4.3.5) 



Pr(A) = (A - Or)pr-l(A) - Pr-2(A), 

with the numbers a,, 6, computed by 



(Pr-l.APr-i) 

■ lIP'-ill 1 



b r = 



llPrll 

IG^r 



r = 1.2 n- 1 . 



(4.3.6) 



(4.3.7) 



This process is numerically stable. 

The only major difficulty encountered by de Boor and Golub lay in the 
computation of the w r eights u^. In seeking to overcome this difficulty, they used 
the reflection of J about its second diagonal. The matrix 



0 0 ... 1 

1 



(4.3.8) 



0 
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is orthogonal and symmetric, so that T 2 = I. It reverses the order of the rows 
and the columns of J, i.e., it transforms J into 

-*n-l <*»-! -K -2 

J = TJT = (4.3.9) 

-hi 
-hi ai 

If, therefore, the elements of J are denoted by a,., 6, then 

Or = "n.l-r. K = 6 n-r- 

The leading principal minors of AI — J are the ti ailing principal minors of AI — J; 
we denote them by p r (A). We prove 

Theorem 4,3.1 For i = 1,2,..., n 

pn-l|Ai)Pn-l(A,) = (blbl ...bn- 1) 2 = b 2 . 

Proof. For once we step out of sequence, and use the notation we will 
introduce in Section 6.2. Let a denote the sequence (2, 3, . . . , w — 1}, then 

P«-i<M = B<aUl), Pn-i(V) = «(<*Un). 

Using Sylvester's theorem (Corollary 2 of Theorem 6.2.2), with B(o) as pivotal 
blodt, we obtain 

0 = B(a) det(B) = 

0 = pn-l(A,)p«_l(Al) - (blt>2...bn-l) 2 ■ 

This result means that the polynomials p l ,(A),p, 1 -i(A),...,pi(A),po(A) are 
the monic polynomials related to the weights 

(4.3.10) 

These weights are more easily constructed than those in (4.3.3). In this proce- 
dure, the terms in the matrix J are computed in the order di,bi,42, . . . ,b„_ i,a», 
i.e., in the order a n ,b n -i ) a n -i....,bi,ai. 

The second method of constructing J is due to Golub and Boley (1977) (133). 
See also de Boor and Saff (1986) (76). It relies on the fact that, once we know 
o(J) and <t(Ji) we may compute the vector xi of first components of the eigen- 
vectors of J; these are the data needed for construction by the Lanczos algorithm 
of Section 4.2. We can carry out the analysis for an arbitrary symmetric matrix 




B(oUl) B(aUl;aUn) 
B(qU l;alln) B(oUn) 
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A € S n , rather than a Jacobi matrix. Barcilon (1978) [19) concentrated on tbe 
eigen vectors corresponding to A* and /i <f rather than using the /i ri to find the 
quantities x ril ; his subsequent analysis did not lend itself to computation. 

If A € then the eigenvalues of Aj are the stationary values of u r Au 



subject to u T u = 1 and the constraint = 0. i.e., u r ej = 
the stat ionary values of 


i 0. Thus they are 


/ = u T Au - Au 7 u — 2i'U T ei t 


(4.3.11) 


where ei = {1,0, ... ,0} and A,i/ are Lagrange parameters. 
/ be stationary yields 


The condition that 


Au — Au — i/©i = 0. 

Since the eigenvectors u, of A span V u , we may wTite 


(4.3.12) 


r 

U = X^o,u,, 

«-l 

and then 

n h 

Au = ^aiAtii = ^ AtQjU., 
*-i «-i 

so that (4.3.12) becomes 


(4.3.13) 



- A)a,u, = i/ei, 

«-i 

and the orthogonality condition u]"u« = Sij gives 

(Aj - A)a,- = imje l = = vx jU 

where we have used (4.2.4). Substituting for a, in (4.3.13) we find 



u 



= uT^r, 

A, — A 



i-1 



(4.3.14) 



and the condition m = 0, and uu = xa, yields the eigenvalue equation 



fM =0 . 



(4.3.15) 



We note that if A is a Jacobi matrix, none of the coefficients Xn will be zero 
(Ex. 3.3.1). The analysis of Section 2.9 shows that the roots (/q)?* 1 of this 
equation will then strictly interlace the (A,)", as in (4.3.1). 

Since xi = {zu,2ai,...,Zni} and xi is the first column of the orthogonal 
matrix X = U T . we have ||x & || 9 = 1 = ^S-itai) 2 , that w have the identity 

^<^)» _nr.-.v-A) 

tr^-A ii;l,(a.-a) 



(4.3.16) 
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(Note that, for large A, both aides approach -1/A.) On multiplying (4.3.16) 
through by (Aj — A) and then putting A = A^ we find 

t = 1,2, ...,n (4.3. 

where 9 indicates that the term j = i has been omitted. The interlacing condition 
ensures that the right hand side of (4.3.17) is strictly positive for each t = 
1, 2, . .. ,n. This equation thus yields x x . 

We stress the importance of the analysis in equations (4.3.11)-(4.3.17). It 
show’s that if A ts an arbitrary symmetric matrix, then <r(A) and <r(A l ) de- 
termine the vector Xi of first components of the normalised eigenvectors of A. 
Conversely, o( A) and x x determine <7(Aj). 

There is a third inverse problem ■which appears in a number of contexts. 
Given two strictly increasing sequences (A,)" and (AJ)" with 

o$a 1 <a;<a 2 <a;<...<a, i <a; ( 4.3.18) 

determine J € 5n such that cr(J) = (At)", and o{3 9 ) = (A?)", where J* = 
(a! — ai)Ei.i -f J. (The matrix J # differs from J only in the 1,1 position.) 

Suppose A € S n is an arbitrary symmetric matrix, and that A* differs from 
A only in the 1,1 position, i.e.. A* = A -f (aj tl — oi,i)Ei.i. We will show that 
<7(A) and <r(A*) determine xi. The eigenvalue equation for A* is 

A*u = Au (4.3.19) 

which we WTite 

Au + K.i - = Au 



<In,a " rfe <*-£>’ 
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where Xn = tijj. The roots of this equation are (A*)J*, *> that 
and therefore 

(“I.i - <*u)4 = (A* - A.) fl ' (£^|) • (4.3.22) 

By comparing the traces of A and A* we see that 

r 

<*!..-a>.i = DA;-Ai)>0. (4.3.23) 

i-i 

Thus, equation (1.3.22) expresses (r*i) 2 in terms of <r(A) and tf(A # ), and the 
interlacing condition (4.3.18) insures that {xn) 7 will be positive. If we know 
that A is a Jacobi matrix then, of course, we can use the Lanczos algorithm to 
determine it. Note that nowhere in the analysis do we need the restriction that 
Ai is non-negative; only the strict interlacing is needed. 

A matrix A is said to be persymmetric if it is symmetric, and also symmetric 
about the second diagonal, the one going from top right to bottom left. Thus 
A is persymmetric if A given by (4.3.9) satisfies 

A = A. (4.3.24) 

If A is tridiagonal and persymmetric, then 

a, = a,,.,.,, b r = 6„_ r . (4.3.25) 

The final inverse problem considered here concerns the reconstruction of a 
persymmetric Jacobian matrix. Now we need only one spectrum, not two. We 
prow 

Theorem 4.3.2 There is a unique persymmetric Jacotn matrix J with a(J) = 
(Ai)*, satisfying 0 $ A, < Aj < • - < A*. 

Proof. The simplest proof is perhaps to show that if the eigenvalues (A,)" 
are known, then it is possible to find the weights for the construction of the 
orthogonal polynomials j>r(A). Indeed if J is persymmetric then the minor 
Pr(A) is equal to JV(A). But then Theorem 4.3.1 shows that 

[p„_i(Ai)| 2 = 6 2 , i.e., p„-i(Ai) = ±b 
so that equation (4.3.10) yields 



wt = ±6/f&(A<). 



(4.3.26) 
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Since the signs of pJ,(A,) will alternate with », then so must the signs in (4.3.26) 
if the w t are to be positive. The magnitude of 6 is irrelevant to the construction 
of the p r (A). See Hochstadt (1979) [182) for another variant of this inverse 
eigenvalue problem. ■ 

Exercises 4.3 

1. Show that if B = A,- 1 - J, then 

B( 1, 2, . . . , n — 1; 2, 3, . . . , n) = (-)*- 1 6 1 4a 

2. Show that the Hi computed from (4.3.22) do satisfy 

•-1 

3. If you like using a computer, then try to reconstruct a Jacobi matrix using 
H aid's method, or that of de Boor and Golub. Start with the matrix J 
with 

o. = 2, bi = 1 , « = 1,2, ... ,n - 1; <u> = 2. 

Set up recurrence relations to give (A,)" and (/r, )i 1 and use these as data 
to reconstruct J. 

4.4 Reconstructing a spring-mass system; by end 
constraint 

We may divide the problem of reconstructing an in-line spring-mass system into 
three stages: 

i) Formulate the problem as an inverse eigenvalue problem for a Jacobi matrix 

J. 

ii) Solve this problem and find J. 

iii) Recover the mass and stiffness matrices M and K from J. 

Stage i) was discussed in Section 3.1; we repeat the analysis here. For an 
in-line system, the frequency equation governing free vibration is 

(K — AM)y = 0. (4.4.1) 

For the system shown in Figure 4.4.1 the matrices K and M are given explicitly 
in (2.2.7). We write M = D 2 , where D = diag(d\,d2 , . . . ,d n ), put Dy = u and 
reduce (4.4.1) to 

(J - AI)u = 0, (4.4.2) 

(4.4.3) 



where 



J = D 'KD *. 
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Stage ii) was the subject of Section 4.3. Given the spectra of the systems 
in Figure 4.4.1a) and b), i.e., <r(J) = (A*)J and <r(J x ) = (/i.JJ* 1 , we construct 
x Xl the \ector of first components of the eigenvectors u, of (4.4.2), and then 
construct J by using the Lanczos algorithm of Section 1.2. 




C ) 

Figure 4.4 . 1 - Two possible wags of constraining the end of a fixed-free system 

It remains to consider Stage iii). By using the explicit form of K in equation 
(2.2.7) we can verify that if 

® = { 1.1 1 ). ( 4 - 4 . 4 ) 

then 

Ke = {*„0,0,...,0}. (4.4.5) 

Physically, this equation states that a static force k t applied to mass mi will 
extend the first spring by unit amount and at the same time displace all the 
remaining masses m2, m3, .... m„ by unit amount to the right, as if everything 
to the right of mi were a rigid body. Since K = DJD we have 



DJDe = DJD) 1, 1, . . . , 1} = {*i,0,0, . . . , 0}, 
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Jd = J {4,4 «U = {fc,/d„o,o,...,o}. 



(4.4.6) 



(Note that D = diag(d x , i^, . . . , d„), while d = («f,, da,-- - 1 d*}.) We need to 
be sure that d so calculated will be a strictly positive vector. We prow 



Theorem 4.4.1 If J € S,, is a non-singular Jacobi matrix, then J -1 is a stiictly 
positive matrix, meaning that each element of J -1 is strictly positive: we unite 
this J" 1 > 0. 



Proof. We use induction. Write 



J = [-b jf]- » = 0}. 

We will haw achieved our goal if we can show that if Jf 1 > 0 then J -1 > 0. 
Suppose 

hi k T 
k H 



J-' = 



then 



so that 



.i-i _ “l ~ bT \ h i kT 1 _ f 1 0 
-b J, J | k H | J 0 1 

— bk T + JiH = I, -bhi + Jik = 0. 



Since J is a non-singular Jacobi matrix, it is positive definite; so therefore is 
J -1 , by Ex. 1.4.2; therefore hi > 0, and so k = Jf'b/u > 0. (Note that the 
product of J," 1 , which is strictly positiw by hypothesis, and the non-negative 
non-zero vector b, is strictly positive.) Therefore, 



H = J I - | bk r + Jf 1 >0. 

(Note that since Jf 1 > 0, all we need in order to prove that H > 0, is that 
bk r > 0, i.e., k > 0; actually though, k > 0.) Thus H > 0. k > 0 and h, > 0 
so that J -1 > 0. ■ 

We may now return to equation (4.4.6). Take the unique reconstructed 
non-singular J and solve 



Jx = J{zi,X 2 ,...,z n } = {1,0,. ..,0} = ei. 



The solution x is strictly positive: x > 0. Thus the solution of equation (4.4.6) 
is 

d = oc, 

for some as yet unknown c > 0. The total mass of the system is 

m = £>,- = £d? = ||d|| J = c^WI 2 . 

»-l »-l 
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Thus, knowing m and ||x|| 2 , we can find c > 0 and d, and thus D. Then 
K = DJD, and because K satisfies (4.4.5), it necessarily (Ex. 4.1.1) has the 
form K = EftE 7 * given in equation (2.2.12), where K = <fu2<?(Ai. A 2 . . .. , A„). 
This completes the reconstruction. 

The reconstruction from the spectra of a) and c) proceeds along similar lines; 
we merely renumber the masses starting from the right (Ex. 1.4.2). 

This reconstruction may be used in a reversed situation: it shows that any 
non-singular Jacobi matrix J may be expressed uniquely as 

J = D l EftE r D- 1 , (4.4.7) 

where D. K are strictly positive diagonal matrices and ! |D|| = 1; this corresponds 
to m = 1 in equation (4.4.6). 

Now' we consider the fixed-fixed case shown in Figure 4.4.2a; there is essen- 
tially only one constraint we can apply, to m,», as shown in Figure 4.4.2b). 
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b) 

Figure 4 4-2 ' A fized-fired system, and a constrained system 



We start our analysis as before. The stillness matrix for the system in a) is 



K = 



Ai -f *2 —A 

-*2 *2 + 



-An 

—An An + A„.fl 



(4.4.8) 



Knowing the spectra (A, )j‘ and (p,), -1 of the systems a) and b) we can construct 
J = D *KD 1 where again M = D 2 Now however 



K{l,l,... ) l} = {* 1 ,0,0,...,0,* n „}: 



(4,1.9) 
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this stales that in order to produce unit static displacements of the masses, we 
must apply two forces, k, at mi and k,,., at m„. Thus 

DJD{1,1,...,1} = kiei +*n« ie„ 

a> that 



Jd = = (* iMJe, + (*n + ,/d„K.. (4.4.10) 

First, consider the equation 

Jy = On- (4.4.11) 

simple algebra shows that the solution is 

Vt =b i b i +i...bn-iP i - 1 /Pi, (4.4.12) 

where Pi is the tth leading principal minor of J (see equation (1.4.6). Since J 
is positive definite, equation (4.4.12) confirms that the solution y is positive, as 
predicted by Theorem 4.4.1. We can find the solution of 

Jx = e, (4.4.13) 



in a similar way (Ex. 4.4.3); all we need here is that, according to Theorem 
4.4.1, x > 0. 

Using x and y we may WTite the solution of (4.4.10) as 

d = (*,/«!, )x + <* B+l /Oy. (4.4.14) 



In particular, 



<*n 



(*i/di)*n + (k m+ i/d n )y H 
kl kn + l Pn-l 

3T* n + 



(4-4.15) 



But P H = nr., and fn-i = 1 17—1* /*,, so that we can write equation (4.4.15) 

<4J16) 

Now consider this equation. The system in Figure 4.4.2a) has 2n-f- 1 parameters. 
Choose one of the parameters, and divide the remaining 2n parameters by it; 
to obtain 2n ratios. The two spectra (AJf and (/i*)?" 1 provide 2n - 1 ratios, 
so one more ratio is needed. The chosen parameter is merely a scaling factor; 
the total mass, or alternatively one individual mass, say m*, would determine 
it. If we take m„ as know*n, then equation (4.4.16) states that the required 2nth 
ratio, fcn+i /m* must be chosen so that 
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This inequality was fiist pointed out by Nyien and tJhlig (1997a) (253). Once 
we have chosen k r +i/m„ satisfying this inequality, then equation (4.1.16) deter- 
mines ki/di, since x n is known, and d n = mj. With k n + \/d n and k\fdi known, 
equation (4.4.14) gives d and hence D and K = D“ l JD . The reconstruction 
is complete. 

The third system is free- free, as shown in Figure 1.4.3a); constraining rrij we 
obtain the fixed-free system in Figure 4.4.3b). 
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b) 

Figure 4*4$ - -4 constraint is applied to a free-free system 

The pair is essentially the same as the pair in Figure 4.4.1, with k\ = 0 . The 
analysis starts as before; the only difference is that the lowest frequency of a) is 
Ai = 0. Still, from <y(J) and <r(Ji) we can construct J uniquely, but now J will 
be singular, i.e., positive semi-definite. 

The stiffness matrix K of system a) will sat isfy' 

K{1, 1, -•-,!} = 0. (4.4.18) 



Now we need a result like Theorem 4.4.1 which covers the case when J is singular. 
It is 



Theorem 4.4.2 If J is a singular Jacobi matrU then the equation Jx = 0 has 
a unique strictly positive solution x satisfying ||x|| = 1 . 

The proof is straightforward; see Ex. 1.4.4. 

Now we may complete the reconstruction. We take J and write K = DJD, 
then (4.4.18) becomes 

Ke = DJDe = DJd = 0. 

Thus d = cx where x is governed by Theorem 4.4.2, and if the total mass m = 1, 
then c = 1. This gives d and hence D and 

K = DJD = EKE r 
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where l< = diag{k l ,k i , . . . ,k n ). Again, we can use this result to show that an 
arbitrary singular Jacobi matrix may be written 

J = D 'EftE'D -1 (4.4.19) 

where now* K has first diagonal entry zero. 

Exercises 4.4 

1. Show that if Ke = k\e 1. and E~* is given by (2.2.10), then KE“ r is 
bidiagonal and E“ 1 KE r is diagonal. 

2. Reconstruct the system of Figure 4.4.1a) from the spectra (At)" and (/x,)? -1 
of a) and c) respectively. 

3. Use the solution (4.1.12) of equation (4.4.11). and the transformation from 
J to J given in (4.3.9) to find the solution to equation (4.4.13). 

4. Provide a constructive proof of Theorem 4.4.2. by writing x in terms of 
the principal minors of J. 

5. Suppose that the eigenvalues (A*)? of the system in Figure 4.1.2a are 

known, as are the eigenvalues (A^)" when the stiffness fc«+i is replaced 
by some unknown stiffness Show that there is a one-parameter 

fcunily of systems, each member of w r hich has the stated eigenvalues. 

6. Show that if J is a non-singular Jacobi matrix, then its inverse C = J 
has the form 

Hit* tilt* ••• til t«n 

tilt* ti2t>2 ••• ti 2 l»n 
tilt>n t i 2 V n ... UnV m 





UiVj X < j 
tljVi i > j 



and that (u g -)y , (t^)J are strictly positive, and satisfy 



— < —<-•• — 
”i *2 i» n ‘ 



This result is quoted in Gantmacher and Kredn (1950) [98]. but may have 
been known earlier. 
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4.5 Reconstruction by using modification 

The simplest way to modify a system is to attach a spring at a free end, thus 
going from the system in Figure 4.5.1a) to that in Figure 4.5.1b). (We have 
renumbered the masses so that the spring is attached at mj.) 




Figure 4.5.1 - .4 spring is added to the system 

This is an example of the analysis of Section 4.3. The spectra for a) and 
b) are <r(J) = (At)" and <r(J # ) = (A*)i respectively. Because we haw added 
stiffness to the system, we have A* > A*, as in (4.3.21). 

i) L T 9e the trace condition to find 






•- i 



ii) L ? 9e 6i = k\/m\ and aj = (k\ + ko)/m\ to find ko/mi = a\ — a 

iii) Use oj — oi and equation (4.3.22) to find x* l9 and hence xi = {m, arai* . , 

iv) Use the Lanczos algorithm to find J. 

v) Use a variant of the analysis given in Section 4.4 to untangle K and M 
from J. 

As an alternative modification we may add mass to the system, specifically 
a mass mj to mi. In this case it is easier to work initially with the original 
equation (4.4.1) than with the reduced equation (4.4.2). Again, w^e start with 
the free-fixed system of Figure 4.5.1a. The eigenvalue problem for a) is 

Ky< = A<My<. 

That for the modified system is 

(4.5.1) 



Ky = AM'y, 
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where M* = M -f mjEjj. Since we have added mass to the system, the 
must 5 



(4.5.2) 



0<AI<A 1 <.-<A,’,<A n . 
Express y as a combination of the y 4 : 

r 

y = £>*• 

M 



Ky = ^o,Ky, = ^A,a,My„ 



(4.5.3) 



i-1 



»-l 



and 



M‘y = 53 aiM y* +*ni E wy» 

i-1 

so that equation (4.5.1) becomes 

n n 

53 A, o, My, = A 53“* My. + AmtE,.,y. 

i-1 i-1 

Premultiply both sides by yT. using the orthonormality condition yTMy ri = 6^: 

ajX, =a,A + AmJy,,i/,, 



and on substituting for in (4.5.3) and equating the first elements of the vectors 
on each side, we find 

" ( w .) 2 



l = Am;£ ^ 



(4.5.4) 



i-1 



In order to use this equation to obtain the first components u u of the eigenvectors 
u, of the reduced equation, for use in the Lanczos equation, we need to express 
y u in terms of ti M . The equation Dy — u gives diy u = u u = Xn that we may 
write (4.5.4) as 



l = XaY a = 

tf A,-A’ m, 

Since the roots of the equation are (A*)" w have 



mj 



trA^A- c n;: 1 (A,-A)- 



(4.5.5) 



Equaling both skies for A = 0 and A — oc, we have 



i-1 



i-1 
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The interlacing condition (4.5.2) gives 

*=n<vA»>i- 

«-i 

The orthonormality condition gives = 1, ao that mJ/nij = a = c— 1 > 

0. Finally, multiplying (4.5.5) throughout by A* — A and then putting A = A g * 
we find 

— AiQX?i = e f[(A; - A*)/n'(A, - A,). (4.5.6) 

i - 1 J-> 

The interlacing condition (4.5.2) ensures that xfi > 0. Now we use xi in the 
Lanezos algorithm, and the untangling procedure as before. 

There are still more ways in which to obtain second spectrum, for which see 
Nylen and Uhlig (1997a) [253|, Nylen and IJhlig (1997b) (251). Ram (1993) (276) 
supposes that the system of Figure 4.5.1 is modified by adding both a mass m 
to mi and a spring ko. He makes use of some simple but powerful results found 
in Ram and Blech (1991) (277). 

We close this section by supposing that an oscillating force F'sinuit is applied 
to the bee end of the spring-mass system of Figure 4.5. la). The matrix equation 
governing the response ysina/t is 



(K-AM)y = Fei. 



Write 



y = ^2<*yi, 

t-1 



where y g is the ith eigenvector, normalised so that 

yT My, = &ij. 



We obtain 
and hence 

s> that 



(A, - A)o, = Fy u , 






(4.5.7) 



When the eigenvalue problem is reduced to standard, J, form, then Dy = u, 
that diyu = ui, = x,i so that we may write 



m _ F *ii 
W mi “ A, — A' 



(4.5.8) 
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where, as usual. Xj = {x 11 , 221 , • . . ,x n \} is the vector of first components of the 
eigenvectors of J. 

The quantity y\/F is called the frequency response function , specifically the 
frequency response function for the displacement yi due to a unit force applied 
at y\. This function may also be identified as a direct receptance for yi, as 
described, for instance, in Bishop and Johnson (19G0) (34). The two spectra 
<y(J) = (A»)J and tf(Ji) = {Pi)\ 1 are the poles and zeros of the response 
function. The interlacing of these two spectra may thus be interpreted as the 
interlacing of the poles and zeros of the response function, a result which is well 
known in control theory. The result of Sect ion 4.3 may thus be stated as follows: 
the response function, and specifically its poles and zeros, uniquely determines 
the matrix J. As we have seen, once we know J and the form of the stiffness 
matrix K. we may untangle M and K from J. See Clad well and Gbadeyan 
(19&5) (106) for an alternative treatment. 

An experimental - theory study of the problem of reconstructing a spring- 
mass system from frequency response data for an actual system may be found 
in Gladwell and Movahhedy (1995) (123) and Movahhedy. Ismail and Gladwell 
(1995) (242). 



4.6 Persy mmetric systems 

It was shown in Section 4.3 that a persymmetric Jacobi matrix J can be recon- 
structed uniquely from its eigenvalues. We shall now consider some physical 
problems relating to persymmetric matrices. Figure 4.6.1 shows a system of 2 n 
masses connected by (2n + 1) springs and fixed at each end. Suppose that the 
system is symmetrical about the mid point , so t hat 

m r =m2«+i-r, hr = Afcn+ 1 -r, (r = 1, 2, . . . , n). (4.6.1) 

The odd numbered principal modes of the system will be symmetrical about the 
mid-point; they will thus be the principal modes of one half (say the left-hand 
half) of the system with the mid-point of the system free, as in Figure 4.6.2(a). 
Thus the odd numbered eigenvalues Ai t Aa,..., A2n-i of the complete system 
will be the eigenvalues of the left-hand half under the conditions fixed-free. i.e. t 

Aai-i = A*, i=l t 2, (4.6.2) 
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Figure i.6.1 - A symmetrical system with 2n -masses 
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On ihe other hand, the even-numbered principal modes of the system will 
be antisymmetrical about the mid-point so that the even-numbered eigenvalues 
A 2 , A t , . . . , A 2f( will be the eigenvalues of the left-hand half under the condition 
fixed-fixed, as in Figure 4.6.2(b). 




(a) 




(b) 

Figure 6.2(a) The odd numbered modes are symmetrical, (b) The even 
numbered ones are antisymmetrical. 



Thus 

A* = A-, 1=1,2 n. (4.6.3) 

This means that the left-hand half, and hence the whole system may be uniquely 
constructed, using the analysis of Section 4.4 from the eigenvalues Ai, . . . , A 2 * 
and the total mass. 

Figure 4.6.3 shows a symmetrical system with 2n — 1 masses and 2 n springs. 
Now the odd-numbered symmetrical modes will be the modes of the left-hand 
half with (m„/2) at the end and free there, as in Figure 4.6.4(a). On the other 
hand, the even-numbered, antisymmetrical modes will be the modes of left-hand 
half with m, fixed as in Figure 4.6.1(b). Thus 

Aai-i=Ai, » = l,2,...,n, 



A 2 .=p„ i= l, 2 ,...,n-L 



(4.6,1) 

(4.6.5) 
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Figure 4.6.3 - .4 symmetrical system with 2n — 1 masses. 




(a) 




(b) 

Figure 4.6.4 - (a) The add numbered modes are symmetrical. 

numbered modes are antisymmetricnl 



(b) The 



4.7 Inverse generalised eigenvalue problems 

In this section we consider how we can reconstruct a finite element model from 
spectral data. 

The eigenvalue problem is 

(K - AM)y = 0, (4.7.1) 

where as in (2.4.10), both K and M are symmetric tridiagonal, K with negative 
codiagonal, M with positive codiagonal. Since one spectrum is insufficient even 
to reconstruct one tridiagonal matrix, it is certainly insufficient to reconstruct 
two. We therefore assume Glad we 11 (1999) (127) that M can be written in terms 
of K: 



M = D J - cK, c > 0, 



(4.7.2) 





4. Inverse Problems for Jacobi Systems 
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where D is an as yet undetermined diagonal matrix with positive entries, and c 
is an arbitrary positive number. Since K has negative codiagonal, M will have 
positive codiagonal. Now 

K - AM = K - A(D 2 - cK) = (1 + cA)K - AD 2 
= (1 +cA)(K — I'D 2 }, v = A/(l + cA). 

Thus (4.7.1) reduces to 

(K - i'D 2 )y = ( 

which, as in Section 3.1, to can reduce to 

(J - v I)u = 0, 

where J = D _1 KD “ l and u = D“*y 

Suppose that (4.7.1) has specified eigenvalues (A*)", where A,* > 0, then J 
has eigenvalues (i'i)i where i /* = A*/(l + cA*) > 0, showing that J, and thus K, 
is posit ive semi-definite. The matrix M can be wTitten 

M = D(I - cJ)D (4.7.5) 

and the matrix I — cJ has eigenvalues 1 — cv i = 1/(1 -f cA, ) >0, showing that 
M is positive definite. 

To reconstruct J we need a second spectrum. If the eigenvalues of (4.7.1) 
under the constraint = 0 are (/i*)"- 1 , then the eigenvalues of (4.7.4) under 
the same constraint will be <r g * = /^/(l + c/ij. We note that the interlacing 

Ai < < Aa < - •• < m„_i < K, (4.7.6) 

yields the interlacing 

i/ 1 <o,<i/ 2 <-.<o n _ 1 <i/„. (4.7.7) 

Having found J. we need to find D so that K = DJD satisfies the characteristic 
stiffness equation (4.4.9). This can be done exactly as in Section 4.4. Gladw^ll 
(1999) (127) finds wider families of systems with the given spectra. See Ram and 
Gladwell (1991) (289) for a different approach to reconstructing a finite element 
model of a rod. 



(4.7.3) 

(4.7.4) 



4.8 Interior point reconstruction 



Suppose, following Gladwell and Willms (1988) (113), we have a spring-mass 
system with n masses, under some end conditions, as in Figure 4.8.1(a). (We 
exclude the free-free condition at this stage.) If a sinusoidal force F sin ujt is 
applied to mass where 0 < m < n — 1, then the response at mass m m * l 

may be calculated as in equation (4.5.7): 






«-l 



A. -A 
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The pole* of this response function are the eigenvalues (Ajy of the whole system, 
A. The zeros of the response function will be the eigenvalues of the system 
constrained so that x mifl =0, i.e., they will be eigenvalues of the systems. B , 
on the left, or C , on the right, of x m +\ % as shown in Figure 4.8.1(b). Different 
ways of assigning the eigenvalues of the constrained system to the two subsystems 
B and C will lead to different reconstructed systems. When this assignment 
has been made, then wo know the eigenvalues (Ajy, (/i*)™ and (i \)\ of systems 
A , B , C respectively; p = n — m — 1. Within themselves these sets of eigenvalues 
must be distinct. There are two cases. 

a) The constrained system has no double eigenvalues. That is, all the 
(/iji'and are distinct; if they are arranged in ascending order and 
relabelled fo)?** 1 , they will satisfy 

Ai <p l < A* <•••<£*_! < A„; 

this is equivalent to the statement that no eigenvector x, of J has a node 
at x m4l , i.e., x m + hi ^ 0 for all i = 1,2, ... ,n. 

b) Two members of a pair ( ) are identical; now* there is an t such that 
Ai = fij = i'k‘y this will occur iff x m +\ f i = 0. There can be more than one 
such pair. 

To analyse the situation we suppose that the eigenvalue equation (4.4.1) has 
been reduced to normal form, (4.4.2). and we partition J as 



ml p 

B ~btn 0 
-bm «m.l -c£,l 
0 -c m ,i C 

where b£ = {0,0,-.-,M, c£., = 




(4.8.1) 
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Figure 4. 8.1 - The mass m m , i is constrained. 









Now we consider the principal minors of AI — J. We denote the leading 
principal minors by p*( A) and the trailing principal minors by The Laplace 
expansions of p„(A) = det(AI — J) using the first m and first m + 1 rows are 

P«(A) = (4.8.2) 

= Pm*i(A)^(A)-£, I p,JA) V _ I (A). (4.8.3) 

We know that 

PW (A) = n(A-A i ), Pm (A)=n(A- Pj ), q f {X) = f[ (A -v k ), 

i-i )-i *— 1 



and thus equation (4.8.2), (4.8.3) give 

Pn(Mi) = — &mPm — l(/* i to,.(/»j) l (4.8.4) 

P-("lr) = (4.8.5) 

In case (a), all the quantities appearing in the latter equations are non-zero, 
so that, apart from the factors 6^ and these equations yield p m -i{Pj) 
4p_i(t/*), respectively. These quant ities are just what is needed to compute the 
matrices B and C, respectively, using Forsythe’s algorithm in Section 3.2. The 
weights for B are given by 



&(«">)* = &Pm-l<Pi)/Pm(Pi). 

= -P"(Pi)/K.(P>)<fr(Pj)|, 

while those for C are 

ftjUi(w*)c = (*'*)/<&(*'*). 

= -Pn(‘ , *)/(^(«'*)Pm(«' l)l- 



(4.8.6) 



(4.8.7) 



To verify that the weights (u'J* are positive, we suppose that p has s v's to 
its left, and p — 8 to its right; then < p . < i /^ M If a number x may be 

written x = (— l) M c, where c > 0. then we say sgn{x) = n. Now w r e can easily 
verify that 



Pm(Pi). %(Pi ) I = fa -j-8,m-j,p-8) 

*) that 

sgn(tuj) t = l + (n-j-a) + (m-j)+p-s 
= 2 n- 2 j- 2a = even 

so that (tt'j)* > 0; we may prove similarly that (tt%) c > 0. 

Thus B may be reconstructed uniquely. At the end, pm-i(A) will be know r n, 
and so the Pm-\{pj) will be knowm. Any one of these values may be substituted 
into (4.8.4) to yield b? ft . The matrix C and i may be found in a similar 



manner. 
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In case (b) there is a common factor A — A, = A — p, = A — in each of the 
equations (4.8.2) and (4.8.3). Cancel this factor and then put A = A,. Since 

Pm+l(A) = ( X ~ <Wl)P-n(A) - 6?„P m -l(A), 

9^.i(A) = (A-a m+1 )q ( ,(A)-^ l9p _ 1 (A), 

we have 

P m .i<P*) = -Om- iGO. 

«p*lO'«r) = -4+1%-iW). 
and thus both equations (4.8.2) and (4.8.3) reduce to 

Pn(^) = (•‘-8-8) 

This is the single equation that replaces the pair of equations (4.8.4) and (4.8.5) 
for the common eigenvalue. Using (4.8.6) and (4.8.7). we may write (4.8.8) as 

= -P^VIpUpX^)] = + C.KIc (4.8.9) 

and we note that \\'\ is a positive quantity. Now we proceed as follow's to find the 
family of Jacobi matrices having the specified eigenvalues. Choose a € (0. j) 
and put 

6 rn (“'» )» = W i a, £(»*)« = IV, sin 2 a. 

If there is more than one triple of common eigenvalues, then this procedure may 
be followed for each. Combine these weights with those corresponding to the 
distinct eigenvalues, and compute B and C. At the final stage pm-i(A) and 
q? i(A) will be known, so that 6& and &£»+ 1 may be found from equations (4.8.4) 
and (4.8.5) for one of the distinct eigenvalues, and there will be at least one, as 
before. 

There is an alternative procedure which elucidates the situation in w r hich 
one or more triples Ai./i,,!/* are equal, and which uses the Lanczos algorithm. 
Express the eigenvalue problem for J in (4.8.1) in terms of the normalised eigen- 
vectors (y^)f 1 and (zj)f of B and C respectively. Thus 

m p 

x = tm(^P,y,) + *m.l<Wl + 9»Z*> 

>-1 k - 1 



where I m = J 1 ^] and = [, 0 ]. The eigenvalue problem becomes 



A-p, -81 




Pi 






Pm 


“«1 •" ~*m A — On,*, “*1 ~tp 




2 m.l 


-i, A-U, 




Qi 


-b x ~»r. 




. 1v . 



(4.8.10) 
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where 





8j = bmVmj, tk = 6m « Ul.A. 


(4.8.11) 


Thus 


(A - M>)P> - = 0. i = l,2,...,m, 

(A-i',)«,-Mm,i=0 I * = l,2,...,p, 




so that 

la . - 


A. + y' + y' w _ o i - 


1 I n 


\“m+l “ 
In v 


7-1 ^ k-l ^ * 

. . . H* ♦ — 1 9 n that 




— / —f t l t t 7~ 

W : V- _ ~P»(A) 

^ 


(4.8.12) 


which yields 


f2 _ 

* K.ivjiqpiPj) h Pm{.»kW p {‘'k) 


(4.8.13) 


for j = = 1,2,...,?, in agreement with (4.8.6), 

may be computed horn 


(4.8.7). Now 






(4.8.14) 



7-1 



k-l 



With (ymj)r ^ (-Ut)i known, from equations (4.8.11)-(4.8.14), B and C may 
be computed by using the Lanczos algorithm. 

In case (b), suppose that there are r > 1 triples (A,*, ('!«}, 9 = 1,2 , ...,r 
such that A,, = /i JV = i-*,. then 






(4.8.15) 



has. as its m+p — r + 1 = n — 2 roots, the n — r non-degenerate A*. In equation 
(4.8.15), ‘means that the degenerate triples are omitted. Now the separate sj 
and tl, and the values W q = sj 9 + tj^ for the degenerate modes will be known. 
Thus as before 



««>, = **', cob 2 a,, f^ = H' g sin 2 a, 

where W q is defined as in (1.8.9). With the a q chosen, the sj q and are ail 
known. Equation (4.8.14) yields and 62,+i as functions of the parameters 
[ct q }i (and note that 6?,. 4*62,+i is invariant) so that the (y«»,») i* and (^i^)f are 
known, and B and C may be calculated from the Lanczos algorithm. 
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An alternative approach to the interior reconstruction problem may be found 
in Nylen and Uhlig (1997a) |253). 

The mass-spring models considered in this chapter are very similar to the 
shear building model used extensively by Takewaki and his coworkers. They 
haw formulated various hybrid inverse problems in which part of a structure 
is given and part is yet to be found in order to yield a structure with specified 
spectral (eigenvalue or modal) properties. Full, definitive description of these 
problems and their use in structural design may be found in the monograph 
Takewaki (2000) (321). Among the original papers most closely related to the 
concerns of this chapter are the following: Takewaki and Nakamura (1995) 

(317), Takewaki, Nakamura and Arita (1996) (318) and Takewaki and Nakamura 
(1997) (319), Takewari (1999) (320). 
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Inverse Problems for Some 
More General Systems 



Words differently arranged have a different meaning, and meanings differently 

arranged have different effects. 

Pascal’s Penstes, 23 



5.1 Introduction: graph theory 

The inverse problems considered in Chapter 4 are special, simply because Jacobi 
matrices are special matrices. In this chapter we will consider &>me slightly more 
general problems but must admit that there are still only a few problems that 
we have been able to solve. 

The special feature of a Jacobi matrix is its structure: it is tridiagonal, 
with strictly negative codiagonal. (It is also positive semi-definite, but that is 
another matter.) The structure of the matrix J in equation (4.1.2) is related 
to the structures of K and M in (4.4.1); K is tridiagonal w'hile M is diagonal. 
The structures of K and M. in turn, derive from the structure of the system, an 
in-line mass system, to which they belong. K, the stiffness matrix, relates to 
the stiffnesses, the connectors, between masses. K is tridiagonal because each 
interior mass m,- 2 < i < n — 1 is connected only to its immediate neighbours 
m^_| and m^j; the end masses m* and each have just one neighbour m 2 or 
m,» respectively. The natural tool for describing and analysing the structure of 
a system is graph theory. 

This is not the place to prove any theorems in graph theory, but it is useful to 
introduce some of the basic concepts. A graph Q is a set of vertices, connected 
by edges l The set of vertices is called the vertex set and is denoted by V; the 
set of edges is called the edge set, £. Figure 5.1.1 shows a graph. This is 
actually an example of a simple, undirected graph. It is simple because there 
is at most one edge connecting any two vertices; the edge connecting vertices t 
and j is denoted by (t f j). The graph is undirected because there is no preferred 
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direction associated with an edge. Henceforth, the terms graph will be used to 
mean a simple , undirected graph. 

The adjacency matrix A of a graph Q is the symmetric matrix defined by 

<k i = 1 i jt j and (», j) € €, 

= 0 otherwise. 



(5.1.1) 



The adjacency matrix for the graph in Figure 5.1.1 is 



A = 



' 0 
1 
1 
0 
0 



1 1 
0 1 
1 0 
0 1 
0 1 



0 0 
0 0 
1 1 
0 1 
1 0 



1 




With any symmetric matrix A we may associate a graph; the rule is 

if i*j then («, j) € £ iff 0. (5.1.2) 

Using this rule we see that the graph associated with a Jacobi matrix is an 
(unbroken) path, as in Figure 5.1.2. The path is clearly one of the simplest 
graphs. 

t 2 5 ••• *n 

Figure 5.1.2 - The graph associated u nth a Jacobi matrix 
Another simple graph is a star on n vertices, shown in Figure 5.1.3. 




Figure 5.1.3 - A star on n t»ertices 
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A (symmetric) bordered diagonal matrix B has a star on n vertices as its 




(5.1.3) 



A periodic Jacobi matrix is one of the form 



ax 6i 
bi aj fci2 





- '• *>n-l 

K — 1 On 



(5.1.4) 



It is tridiagonal except for the terms bn in the top right and bottom left. The 
underlying matrix is a ring on n vertices as shown in Figure 5.1.4. 




Figure 5.1.4 - -4 ring on n vertices 



The graph 
tion 2.3 in the 
5.1.5. 



with a pentadiagonal matrix, such as occurred in Sec* 
of the vibration of a beam, is a strut, as shown in Figure 




Figure 5.1.5 - The strut on n (even) vertices is the underlying graph of a 

pentadiagonal matrix 
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The graph associated with a 2 x 2 block tridiagonal matrix is also a strut, 
but now one with double connections, as shown in Figure 5.1.6. 




The graphs shown in Figs. 5. 1.1-5. 1.6 are all connected graphs: there is a 
chain consisting of a sequence of edges connecting any one vertex to any other 
vertex. Note that the intersections of the diagonals in Figure 5.1.6 are not 
vertices of the graph. 

The graphs shown in Figure 5.1.7a), b) are disconnected. 




b) 



Figure 5.1.7 - Renumbering does not essentially change a graph 



In order to test whether the underlying graph of a given (symmetric) matrix 
is connected or not, we note that renumbering the vertices of a graph does not 
change the essential character of a graph; the graphs a) and b) in Figure 5.1.7 are 
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essentially the same. Renumbering the vertices of a graph leads to a rearranging 
of the rows and of the columns of any (symmetric) matrix based on that graph. 
When a graph is disconnected, it may be partitioned, as in Figure 5.1.7a) into a 
set of connected subgraphs. Then we can always rearrange the numbering, as 
in b) so that vertex numbers in any one connected subgraph form a consecutive 
sequence. The adjacency matrices of the graphs a) and b) are 



0 


0 


0 


1 


0 ' 




' 0 


i 


0 


0 


0 ' 


0 


0 


1 


0 


1 




1 


0 


0 


0 


0 


0 


1 


0 


0 


1 


. A 2 = 


nr 


0 


TT 




— I 


1 


0 


0 


0 


0 




0 


0 


1 


0 


1 


0 


1 


1 


0 


0 

• 




0 

. 


0 


1 


1 


0 



We see, in this example, that when the vertices are renumbered so that each 
connected subgraph has consecutive numbering, then the adjacency matrix splits 
into two separate submatrices: such a (symmetric) matrix is said to be reducible. 
A symmetric matrix A is said to be irreducible iff it cannot be transformed to 
the form 

<5.1.5) 

by any rearrangement of rows and columns. If it is reducible, then it can be 
transformed to the form (5.1.5), and of course B and C may perhaps themselves 
be reduced further. Note: The concepts of connectedness of a directed graph, 
and the corresponding concept of irreducibility of a general (not necessarily 
symmetric) matrix, are more complex than those described here. See Horn and 
Johnson (1985) [183] Section 6.2.21. 

Now we may state the general result. 



A = 



B 


LJL 




1 



Theorem 5.1.1 The (symmetric) matrix A is irreducible iff its underlying graph 
is connected. 



It is easy to check that if a spring (other than *i) is removed from a spring 
mass system such as that in Figure 4.4.1, then the underlying graph becomes 
disconnected, and the stiffness matrix becomes reducible. 

A tree is a special kind of connected graph: one which has no circtitfs. Now 
there is a unique chain of edges connecting any one vertex to any other. The 
path and the star are both trees, but a ring, see Figure 5.1.4, is not a tree. A 
connected graph has one or more spanning trees. If Q is a connected graph with 
vertex set V, then a spanning tree S of Q is a maximal tree with the vertex set 
V; if any more edges in £ were added to «S then it would cease to be a tree: it 
would have a circuit. Figure 5.1.8 shows three possible spanning trees for the 
graph Q in Figure 5.1.1. 
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It may be proved that all the spanning trees of a given graph Q have the 
same number of edges. 

Kabben (2001) [2-13], in a wide ranging paper, discusses Green's matrices for 
trees. 

5.2 Matrix transformations 

In the first part of this book we are concerned very largely with matrix eigenvalue 
problems. One of the basic questions we face is this: ‘What operations. i.e. f 
transformations, may we apply to a matrix, or a matrix pair, which will leave 
its eigenvalues unchanged, i.e., invariant? 1 We now discuss this question. 

Suppose C, A € M n . The set of matrices C - AA is called the matrix pencil 
based on the pair (C, A). As stated in Section 1.4. the eigenvalues of the pair 
(C. A) are the values of A for w f hich the equation 

(C - AA)x = 0 

has a non-trivial solution x € V u . The eigenvalues are the roots of 

det(C - AA) = 0. 

Suppose P. R € M n are constant matrices, i.e., they are independent of A. Since 

det(PCR - APAR) = det(P) • det(C - AA) • det(R) 

we may deduce that if P. R are non-singular, so that det(P) / 0, det(R) ^ 0, 
then 

det(PCR - APAR) = 0 iff det(C - AA) = 0, 

so that the transformation ‘premultiply by P, and post multiply by R’ leaves the 
eigenvalues invariant. The transformation is called an equivalence transforma- 
tion. It is a special equivalence relation (Ex. 5.2.1). 

In general, an equivalence (transformation) will transform a symmetric pencil 
into an unsymmetric pencil. Those which preserve symmetry are characterised 
by 



P = R r . 



(5.2.1) 
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An equivalence changes a pencil A — AI into PAR - APR If it is to change 
A - AI into B - AI, then we must choose P, R so that PR = I, i.e., 

P = R~\ (5.2.2) 

An equivalence with this property is called a similarity (transformation). An 
equivalence which satisfies both (5.2.1) and (5.2.2) is called a rotation or an 
orthogonal transformation. We reserve the symbol Q to denote the l P‘ of such 
a transformation. Equations (5.2.1), (5.2.2.) show that 

QQ r = Q r Q = I : (5.2.3) 

Q is an orthogonal matrix ; the matrices U and X in Section 4.2 were orthogonal 
matrices. We recall that the columns (row’s) of an orthogonal matrix are mu- 
tually orthogonal, and each column (row f ) has norm 1; if Q = 
then 

q T<b = V (52.4) 

If n = 2. an orthogonal matrix has the form 

*-[Zl -*/]■ <“*» 

When n = 2, the eigenvalue problem relates to a plane, and this Q corresponds 
to a rotation of the x, y axes through an angle 0 about the r-axis. 

It is difficult to write down the most general expression for an orthogonal 
matrix in Af„. Instead, we use the fact that a product of orthogonal matrices 
is itself orthogonal (Ex. 5.2.3). 

There is a particularly simple and powerful orthogonal matrix which can be 
constructed by making a rank-one change to the identity matrix: 

Q = I - 2/ixx r (5.2.6) 

will be orthogonal if 

QQ r = (I - 2pxx r )(I - 2pxx T ) 

= I - 4pxx r + V(x r x)(xx r ) = I, 



i.e., if p is chosen so that 

p = l/x T x. (5.2.7) 

Such a transform is called a Householder Iransjotmalion ; note that Q in (5.2.6) 
is symmetric, i.e., Q = Q T . 

Householder transformations are used in various contexts; one is the reduc- 
tion of a symmetric matrix to tridiagonal form, as we now describe. 

Suppose Q is given by (5.2.6), and A € S„. We wish to choc** Q, i.e., find x, 
so that the transformed matrix QAQ = B has zero elements in its first row and 
column, except for the first two, #>11,612. First consider the postmultiplkation 
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by Q. and use the abbreviation p,, for row 1 of a matrix. With Q given by 
(5.2.6), we haw 

C = AQ = A - 2p(Ax)x T 

p,(C) = p,(AQ) = Pi<A) — 2p(p,(A)x)x T 



Thus c u = a u - 2p(p,(A)x)x il i = 1,2 n. 




We now choose 




x i = a u , « = 3, 4, . . . , n. 


(5.2.8) 


Then a,- =0, i = 3,4, . . . ,n if 




2p(p,(A)x) = 1. 


(5.2.9) 


This gives one equation for the remaining unknowns xi, xj. 
premultiplication: 


Now carry out the 



so that 

p I (B) = p 1 (C)-2pp l (x)(x T C). 

Thus if the premult iplicat ion is not to change the zero elements in the first row 
of C, we must chocs* x x = 0. Now equations (5.2.7)-(5.2.9), give 

2 (u 12 j 2 + a? 3 + ■ • • + a 2 ln ) = x* + (n? 3 + ■ • • + o?„), 



which yields 


i 2 = a l2 ± S, 


(5.2.10) 


where 


■ ■ - 

n 






(5.2.11) 




t-2 




Thus the required x is 


x = (0, a l2 ±S, a l3 ,...,a ln } 


(5.2.12) 



and for numerical purposes we choose the sign of 5 to be that of o 12 . 

This is the basic Householder transformation; it reduces an arbitrary sym- 
metric A to a matrix 

B =[b“ £]• <“•“> 

where b = &iei. This completes the first step in the reduction to tridiagonal 
form. Now we apply another Householder transformation to the submatrix 
Bi, using a new x with x\ = 0 = xa. This second transformation will leave 
an, b, b r unchanged, and will eliminate all but the first tw r o elements of the first 
row and column of Bi . After n — 2 applications, the matrix becomes tridiagonal. 
Once the matrix has been reduced to tridiagonal form, its eigenvalues can easily 
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be located by using the sign count function s r (A) of Section 3.1. Details on 
the numerical implementation of this reduction may be found in Bishop, Glad* 
well and Michaelson (1965) (33] (Chapter 9), Golub and Van Loan (1983) (135] 
Section 8.2. 

We make two comments. Because z x = 0, the Q in (5.2.6) may be written 




(5.2.14) 



where Qj is an orthogonal matrix in This has an important consequence. 

Not only does the transformation preserve <r(A), i.e., o( A) = <r(B), but also 
<7(Ai) = tf(Bi). 

Secondly, we can use a trivial modification of the Householder transformation 
to reduce a general symmetric matrix A to, say, pentadi agonal, form. We take 



x = {0, 0, a 13 ± S, a Mi . . . , a ln } (5.2. 15) 

where S 2 = This transformation preserves <r(A) f <7(A 1 ),<r(A lt 2), 

where the last denotes the spectrum of A with rows and columns 1 and 2 re- 
moved. 



1. An equivalence relation , ‘o is related to b\ written o/?6, has three defining 
properties: 

• reflesivity . qRq 

• symmetry . if aRb then bRa 

• transitivity , if aRb and bRc. then qRc 

A set of elements related by an equivalence relation is called an equhalence 
class. Use the joint operation ‘premultiply by P and postmultiply by R' (with 
P, R non-singular) to define an equivalence relation and an equivalence class for 
matrix pairs (C, A). 

2. Show that the transformation B = QAQ r defines an equivalence relation 
and a corresponding equivalence class. 

3. Show that if Q 1 .Q 2 are orthogonal, then so is Q 1 Q 2 . Show by coun- 
terexample that if Qi,Q 2 are symmetric, then Q 1 Q 2 is not necessarily 
so. 

4. Show that if x is given by (5.2.11), then p in (5.2.7) is given by 25(5 4- 
tfi2)/i = 1- 

5. Verify that the Q obtained as a result of n — 2 successive Householder 
transformations has the form (5.2.14). 
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5.3 The star and the path 

In Section 5.1 we noted that the graph associated with a bordered diagonal 
matrix (5.1.3) is a star on n vertices, as in Figure 5.1.3. There is a particularly 
simple inverse eigenvalue problem connected with a bordered diagonal matrix 
B: construct B so that <r(B) = (A, <r(Bi) = (p,-)" -1 - The usual variational 
arguments show that the two spectra must interlace, at least in a loose sense: 

A!<P!<A 2< • (5.3.1) 

For simplicity we assume that the are distinct: 



Pi < Pa < • * • < Pn-l- 


(5.3.2) 


We write B in the form (5.1.3), i.e., 




-[? m]' 


(5.3.3) 


where M is diagonal, and f> = Clearly, we can make <t(Bj) = 

(/0?” 1 by taking M = ^ 2 ,..., p n -i)- The trace condition giv^es 


n w-l 

»-i «-i 


(5.3.4) 


Now consider the eigenvector equations for B: 




6,1>1 + (p, - A)» i+1 =0, > = 1,2 n - 

(ai - A)m + kv,*i = 0, 


- 1, 


which give the eigenvalue equat ion 




t%-l jo 




This is to haw roots (A,)?, so that 




. . ^ S nr.,(A-Ao 


(5.3.5) 


and hence 

,2 ~ nj'-il/*. ~ A j) 

* = , " u 


(5.3.6) 



where, as usual. 1 denotes i / j\ the interlacing condition (5.3.1) yields 5 * > 0. 
We can choose the sign of 6j to be + or Because we have assumed that 
the /i, are distinct, a given /i, can coincide only with its neighbours A i or A<+i. 
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Equation (5.3.6) shows that b 4 = 0 iff /i- coincides with either of these two A ? s. 
If = 0, then the edge (l,s + 1) is absent from the underlying graph. 

Having constructed the bordered diagonal matrix B. we have a new way to 
construct a tridiagonal J such that <r(J) = (Ajy, <7(Jj) = (/i,-)"” 1 : we can 
apply Householder transformations to B to get J. On account of Ex. 5.2.5, the 
transformation will have the form 




On carry ing out the multiplication, we find 



QT j iQi = M. QlV, = 6. 



(5.3.7) 



(5.3.8) 



(5.3.9) 



The first equation shows that the eigenvectors of Ji are the columns of Qi: the 
tth eigenvector is 

q. = {«<.«» «(.-!»}• (5.3.10) 

The second equation shows that, apart from the factor b ly the vector 6 is the 
vector of first components of the eigenvectors of J i : 

& = (5.3.11) 



Thus, apart from the factor 6 lf & is the vector Xj needed for the construction 
of Ji from the Lanczos algorithm of Section 4.2. The factor is given by 

= ii&ii- 

Note the difference between (5.3.6) and (4.3.17): the former, according to 
(5.3.11), gives the first components of the eigenvectors of Jj; the latter gives the 
first components of the eigenvectors of J. 

Sussman-Fort (1982) [312] discusses connections between the inverse eigen- 
value problems for Jacobi and bordered matrices. 



Exercises 5.3 

1. Explore what happens to J w r hen one or more of the A ? s coincides with a 
P- 



5.4 Periodic Jacobi matrices 

In Section 5.1 we showed that the graph underlying a periodic Jacobi matrix is 
a ring on n veitices.The following analysis is due to Ferguson (1980) (87) and 
Boley and Golub (1984) (35), Boley and Golub (1987) (36). 
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A periodic Jacobi matrix J pcr has 2 n terms, We show how to 

construct J^ r from <r(J f * r ) = (**)?» <’(j pir .i) = (/i,*)"” 1 and one extra piece of 
data; 

0 = tnb2...t>n. (5.4.1) 

It is convenient to consider two matrices, the original matrix J yer of (5.1.4), and 
another matrix J^. with 6,» replaced by — bn- We suppose <r(J^ r ) = (A, r )" ; 
dearly there are relations betw r een the X i and the A,. The A* and will again 
interlace as in (5.3.1), as will the Ar and p it again w r e suppose that the (p 4 )|~ l 
are distinct, i.e., (5.3.2) holds. 

We start by constructing two bordered diagonal matrices. B horn (A*)y and 
(Wj)*" 1 , B“ from (\“)y and (m*)*" 1 . They will have the form 



B = 


a, b r ‘ 


. B" = 


\P, 

i 

H 




b M 


> 


b* M 1 



(5,1.2) 



Here aj,S will be given by (5.3,1), (5.3.6), and of, 6“ will be obtained from 
(5.3.4), (5.3.6) by replacing A, by Af. 

Since «r(Jp„.) = o(B) and <x( Jpcr.i) = cr(B,), J„ r and B are related by an 
orthogonal transformation of the form 



B = 



<*i 

6 



&r l = 

M J 



0 

QT 



a, Mf + 

6,ei + 6„e n _i A! 



1 

0 



0 

Q. 



where ©i = {1,0 ...., 0}, e n -i = {0,0 ,...,!} are in Vn-i and similarly 



[•»“ 6 ' T 1 = [ 1 ° 


of -6„eX_, 


1 6- M | l o Qf 


Mi ~ A, 




The subblocks of these equations corresponding to b and b are 



6 = Qf(6,ei + &«e„_i) 

6- = Qi'(b|©l - &n e w-l) 



which on addition and subtraction give 



6 + 6- = 26,Qfe, = 26,x, 

6-6- = 26„Qfe n _ 1 = 26,x n _ I 

where, as in (5.3.9), xi b the first column of Qf and x„_i is the (n - l)th 
column. If w r e know* b and b“, then these equations give b\ and bn (up to sign) 
since ||xi|| = 1 = ||x„-i||. Once w have ai,&i and xi we may compute Jp«r.i 
from the Lanczos algorithm as before. 

However, in finding B and B^, specifically in finding b and b”, we assumed 
that we knew both the (A «)? and the (Af)i- We complete the analysis by 
showing that we can in fact find b~ from the (A*)" and in (5.4.1). 
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The periodic Jacobi matrix J^ r differs from a regular Jacobi matrix only in 
the presence of the entries b r in the corners. This means that det(AI — J^ r ) and 
det(AI — J^.) will differ from the nth principal minor p n (A) only by quadratic 
terms in In fact 

del(AI — Jpcr) = ni,(A-A0 = p H (A)-6ir„_ 2 (A)-2fl 
detlAI-J^,) = nr-i(A-Ar) = p„(A) - 6*r n _ a (A) + 20 

where r n -2 is the principal minor taken from row’s and columns 2, 3, . . . ,n — 1. 
Subtacting these two equations, we find 

n<A-AT) = n(A-A,) + «l/i 

•-1 *-l 



This means that we can express (6 ; - ) 2 in terms of (A,)” and (/^)V 



n-1 



, ; _ ,2 _ n"-i 0 *i - K ) _ II?-1 + 

i ir_v '<*-*>“" nr.v ^ - *> ' 

But this expression is not automatically non-negative if the (A*)" and (/i,)" 1 
satisfy the interlacing condition. We must examine this more closely. Suppose 
first that 0 = 0. The expression is certainly non-negative, and actually positive 
if the A ? s and p’s strictly interlace. If they strictly interlace then, from continuity 
considerations we can conclude that, for each value of j, the expression will be 
non-negative for 0 lying in a closed interval (— ej> fj\ around zero, e ^ > 0, fj > 0. 
This means that all the (6j) 2 will actually be non-negative in the intersection 
of these closed intervals. For 0 in this intersection the problem as posed, has a 
dilution; for 0 outside this interval it has no (real) solution. Boley and Golub 
(1987) [36) present an algorithm to compute in this way. See also Boley and 
Golub (1984) (35). Xu (1998) [339) provides a detailed analysis of the problem 
and shows (Theorem 2.8.3) that there is a solution iff 



it, \*h ” **l £ 20(1 + (-1)*” >+I ) <5.4.3) 

*-l 

for all j = 1,2,..., n - 1. Note that if (A,*)y and (/i # )"“ l are given, then the 
inequality (5.4.3) provides an upper bound for 0. Andrea and Beny (1992) (9) 
provide a completely different approach to the problem via continued fractions. 



5.5 The block Lanczos algorithm 

In Section 5.1, we exhibited Figs. 5.1.5 and 5.1.6, and showed that the matrices 
underlying these graphs were pentadiagonal or block tridiagonal. In order to 
develop methods for living inverse problems for such systems, we need a block 
version of the fundamental Lanczos algorithm described in Section 4.2. 
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First we recall the original scalar version: Given a symmetric matrix A, and 
a vector X[, compute a Jacobi matrix J as in equation (4.2.1) and an orthogonal 
matrix X = (x ll x 2 ,...,x n | such that A = XJX T . 

The algorithm proceeds by using the two equations 

J = X r AX. AX = XJ, (5.5.1) 

alternately. Thus the (1.1) term in (5.5.1a) gives 

a, = xf Ax, 

and the first column of (5.5.1b) gives 




and so on. 

We now construct a block version of these equations, following 
Boley and Golub (1987) (36] . We start with a symmetric matrix A € S’„ and 
suppose n = pa for some integer a. We will reduce A to a block tridiagonal 
matrix J, where 

A, -Bf 
-B, A 2 -Bl 

-B._, A, 

Here Ai,...,Aj are symmetric, i.e. t in S Pi and the B, are upper triangular 
matrices in A f p . We assume that in addition to A, we are given p oithonormal 
vectors (x,)j € V n which form the columns of Xi = [xi,X2»... , x p ) € j 
T he matrix Xi therefore satisfies XfXi = I p . 

The aim of the procedure is to construct J and an orth<>gonal matrix X = 
(Xi t X 2 ,...,X.) such that 

A = XJX r . 
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Just as in the scalar Lanczos process, we consider the two equations 

J = X r AX, AX = XJ. ( 5 . 5 . 4 ) 

The first pxp block of the first equation gives 



while the first n x p block of the second gives 



AXi = XiAi — XjBi 



which we rewrite as 

X 2 B! = Xi Ai - AXi = Zc- 

In the scalar version we had 61X2 = Zj, from which we immediately concluded 
that 61 = HzjII, and hence x 2 = In the block version we have constructed 

Z 2 € My and we wish to write it as X 2 B,. Write X 2 = |y,, y 2 , . . . ,y p ], Z 2 = 
••>**] and 

[ 611 612 ... 6i p 1 




then finding (y«)i and the elements of Bi is essentially a Gram-Schmidt process: 
finding orthonormal combinations of the vectors (z*)f. Thus 

friiyi = 21 implies bn = ±||zi|| f yi = z\/b\\ ( 5 . 5 . 5 ) 

and then 

fri2yi 4 - b&y2 = z 2 

gives 

612 = yfz 2 , 622 y2 = Z2 - 6i2yi = w 2 

so that 

fr22 = ±||w2|| t y 2 = W2/622 etc. ( 5 . 5 . 6 ) 

The Gram-Schmidt process is closely related to the QR algorithm. The de- 
composition X 2 Bx = Z 2 involves writing Z 2 as the product of X 2 which is in 
but which satisfies XjX 2 = Ip, and an upper triangular matrix Bj € Af p . 
Because X 2 is not simply an orthogonal matrix in A f pi the usual QR algorithm 
has to be modified to effect the decomposition. 

Now we can proceed as before. We have found X 2l so that 



and 



AX 2 = — XiBf + X 2 A 2 - X,B 2 
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so that 



X 3 B, = -X.Bf + X 2 A, - AXj = Z, 



from which X3, B 2 may be found, as before, by Gram-Schmidt. Note that 
different choices for the square roots, as in (5.5.5) and (5.5.6) will lead to different 
matrices J. Boley and Golub (1987) (36) present a detailed algorithm for the 
process. 

Further studies on the block- Lancz os algorithm have been carried out by 
Underwood (1975) (325) and Golub and Underwood (1977) (134). See also 
Mattis and Hochstadt (1981) (222). A completely different and highly efficient 
procedure for the solution of band matrix inverse problems has been developed 
by Biegler-KOnig (1980) (28), Biegler-Konig (1981a) (29), Biegler-KOnig (1981b) 
(30), Biegler-KOnig (1981c) (31). See at*> Gragg and Harrod (198-1) (153) for 
a procedure based on Rutishauser's algorithm; they explore the connections to 
a number of other problems. See also Glad well and Willms (1989) (114) and 
Fiiedland (1977) (92), FMedland (1979) (93), and particularly, Ghu (1998) (58). 



5.6 Inverse problems for pentadiagonal matrices 



We could pose an inverse eigenvalue problem for a general symmetric matrix 
with 2p + 1 bands, as in Boley and Golub (1987) (36). Instead, we will confine 
ourselves to the case p = 2, a pentadiagonal matrix A. The pentadiagonal case 
occurs in the inverse problem for a vibrating beam, but we shall defer considering 
the beam until we have discussed positivity in Chapter 6; the pentadiagonal 
matrix giving the stiffness matrix of the beam has a very special form; certain 
terms in it must be positive, and others must be negative. In this section we 
will not be concerned with these matters of sign. 

Suppose we are given 

«r(A) = (A.)™, <r(A,) = (p.)?- 1 , o(A, <a ) = W 3 , (5.6.1) 

where, as before, <7(A 12 ) denotes the spectrum of A 12 when its first two rows 
and columns are removed. Clearly the eigenvalues must interlace; and for 
simplicity we assume that the interlacing is strict. 



Ai < < Aa < a * a </*,,-! <A„, (5.6.2) 

Pi < 1/1 <p 2 < ••• < |/*- 2 < Mn-l- (5.6.3) 

Our aim is to construct A such that (5.6.1) holds. We WTite 




b T 

Ai 



(5.6,1) 



where only the first two components of the vector b are non-2ero. We denote 
the eigenvector matrix of A by Q r and of A x by Q<*> so that 



Q r AQ = A, Q^AiQW = M. (5.6.5) 
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The eigenvectors of A are therefore q,, where Q = [qi,q2* • •• ,q»] while those 
of A, arc qW where Q<» = |q^ qi‘2,]. 

We start by constructing a bordered diagonal matrix, as in Section 5.3: 




(5.6.6) 



such that <r(B) = (A t )J, and *(M) = The term is given by the 

trace: 

a i = < 5 - 6 - 7 ) 

t-i «-i 

while 6 is given by (5.3.6): 

Now, following equation (5.3.8) we relate A to 



r o, b T ■ 




1 0 


ai b r 


1 0 


[ b M 




0 Q<'> T 


b A! 


0 Q(» 



(5.6.9) 



As in (5.3.9), we have 



Q"» T b = 6 



(5.6.10) 



Now however, in eontjast to the situation in Section 5.3, b is not just a multiple 
of ej, so that b does not give the vector of fust components of the eigenvectors 
of A( . But we can use the analysis of Section 4.3 to obtain the first components 
of the eigenvectors of A and A t : 



, - a.) 

9u ■ n; 1 -, 



ri’-i ("> - 

,91 ' ) " H7-; *(mj _ Mt)' 



(5.6.11) 



To apply the block Lanczos algorithm to construct A we need not just the 
vector xi of first components of eigenvectors of A, but also xa of second com- 
ponents, making up Xi = [X 1 .X 2 ] € Partition the vector q*: 



*-[n 



y. € v._,. 



(5.6.12) 



Since q, is the ith eigenvector of A, and A is given by (5.6.4), we may write 



<j| b T 

b A, 




(5.6.13) 



so that 



9iib + Aiy, = A«yi. 




no 
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Now premultiply by Q ,l,r to obtain 

9i,Q |,|r b + Qi l,T Aiy, = A,Q"» T y„ (5.6.14) 

But equation (5.6.10) gives Q a ' T b = b, and equation (5.6.5b) gives Qi" T Ai = 
MQ (1,T , so that equation (5.6.14) gives 

9ll 6 = -(M-V)Q (,)T y* 

and hence 

y* = -fli*Q ,l >(M - A.ir'b. 

We need just the first term in y,*; it is 

w-i Ji); 

yn = ' ' J [ '• ‘=1,2 (5.6.15) 

i- 1 Pi ~ 

Since Sj is given by (5.6.8), and are given by (5.6.11), this equation yields 

p„, and hence x 2 = {p n ,ifi 2 ,. . - , !/,„). 



Exercises 5.6 

1. Verify that the vector x 2 given by (5.6.15) is indeed orthogonal to x>, as 
required. 

2. Extend the procedure described in this section to the general case of a 
2p + 1 band matrix. 



5.7 Inverse eigenvalue problems for a tree 

The inverse eigenvalue problems for a path and a star are particular examples of a 
general problem. Both the path, as shown in Figure 5.1.2, and the star, in Figure 
5.1.3, are trees, as defined in Section 5.1. The matrices corresponding to these 
trees are a Jacobi matrix J, or, as we will choose here, a sign-reversed Jacobi 
matrix A = j, and a bordered diagonal matrix respectively. In both problems, 
two spectra were specified, namely <r(A) = (A,-)?, and <r(A,i ) = (/q)? -1 ; the 
second spectrum corresponded to the eigenvectors u set to zero at a prescribed 
vertex, vertex 1. In both cases the two spectra had to satisfy' the Cauchy 
interlacing inequalities 



A. <#i l <A 2 <...< / q,_,<A n . 



(5.7.1) 



In both cases also, if the inequalities (5.7.1) were all strict, the matrix A was 
ineducible, and the corresponding graph Q was connected. 

The purpose of this section is to sene as an introduction to an important 
paper by Duarte (1989) [81]. This paper reviews the history of inverse eigenvalue 
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problems for trees, and establishes a general result. We will present analysis 
cowring the simpler parts of the general case. As we will do sometimes in 
Chapter 6. Duarte labels eigenvalues in decreasing order, and we do the same. 
Specifically, we will show' that if Q is a tree on n vertices V, and if two spectra 

(\)UHi)i~ l “« g> wn - satisfying 

Ai > > X 2 ■ ■ ■ > /i„_i > An > 0, (5.7.2) 

then we can find a symmetric matrix A € S n on G such that o = <r(A) = 
(AJf, <*i = ^(A ?1 ) = (p,)?" 1 . We take the strict interlacing and the positiv- 
ity condition for simplicity; Duarte relaxes these conditions. 

We start by observing that the two cases that we have considered so far, the 
path (Jacobi), and star (bordered diagonal), have common features. First, we 
note that the entries of the constructed matrices may be considered as functions 
of the data a*0\. Secondly, we note that in both matrices there are n 2 — n — 
2 (n — 1) = n 2 — 3n -f 2 constant functions, w'hich in fact are all zero. This 
suggests the following questions: 

1. Can the constant functions appearing in A be other than the zero function? 

2. Can the number of these constant functions be increased? 

The answer to the first question is NO. For if A € S n has eigenvalues (A ,*)J 
with maximum modulus A, then (Ex. 5.7.1) |atj| < A, so that A can have no 
fixed entry, independent of the eigenvalues, other than zero. To answer the 
second question w r e note that if the inequalities (5.7.2) hold, then A must be 
irreducible. For if A were reducible, i.e., after possibly renumbering the vertices, 
it could be written 




then A and A.j (w'hich after renumbering, would be A n ) would have a common 
eigenvalue, a situation that is precluded by (5.7.2). Thus A is irreducible and G 
is connected. Now* we note that A must be positive definite so that no diagonal 
term a ii can be zero. The maximum number of zero entries will be attained for 
matrices whose graph is a tree, and this number is precisely n 2 — 3n 4* 2 (Ex. 
5.7.2). Thus the answer to the second question is NO also. 

Having answered these questions, we proceed to the analysis. We start by 
considering a tree G, choose a vertex of V, label it 1, and see the effect of deleting 
vertex 1 - this is the graph corresponding to deleting row 1 and column 1 of A. 
First, we need a symbol, A/", to denote the set of m vertices j of Q which are 
connected to vertex 1. Now' we use G' to denote the graph obtained from Q by 
deleting vertex 1. Figure 5.7.1 shows two examples. In Figure 5.7.1a, where 
vertex 1 is at the end of a path, j\f = {2} and G* is the connected graph with 
vertices {2,3, 4,5,6}. In Figure 5.7.1b, A r = {2,4} and G' has two connected 
components, one on either side of vertex 1 ; we call these G 2 » Ci respectively. In 
general, Q f will have m connected components which we label Gj , j € the 
corresponding matrix A,i will have m irreducible components. 
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i — a — > — j — a — i i i i i i i 



i l i i J i—i i—i—i 



a) b) 

Figure 5.7.1 - Deleting vertex 1 from a path. 

Figure 5.7.2 shows another example, a star. Now M = {2. 3, 4, 5} and Q' has 
4 connected components: SJ = {j}, j = 2, 3, 4, 5. 



•5 



• 

• 3 

2 3 2 

Figure. 5.7.2 - Deleting the centre of a star. 

Finally, we need a symbol for ihe graph obtained by deleting vertex j 6 A r 
from </j; we call it S”. Figure 5.7.3 shows these subgraphs for the graphs Q' in 
Figure 5.7.1. 




• • • • • • • 

3 4 5 6 3 5 6 

a ) b) 

Figure 5.7.3 - The subgraphs Q" for the graphs & . 

Note that for the star, the vertex set of (% is empty because the vertex set 
of £ is 0). 

Having established notation that allows us to see what happens when we 
delete a vertex of a graph, we need to consider the two sets of eigenvalues, and 
how these relate to the matrix A. To do this we first return to two examples we 




5. Inverse Problems for Some More General Systems 



113 



haw already Healed, those corresponding to deleting an end vertex of a path, 
and the centre of a star. 

First, the path with end vertex 1. The eigenvalues (A,)? and (p,)J -1 are the 
zeros of the trailing monic principal minors, /’'(A), J* ,(A) respectively, and, in 
the notation of equation (4.3.4), these are linked by 

KM = (A - a,)^_,(A) - t$K-2( A). <5.7.3) 

We note that the graphs corresponding to /?„ K-i< K - 2 are precisely Q, Q' = Q' 2 
and Q" = Qji in fact K< K-i, K - a are the characteristic polynomials A of A, 
and of the submatrices of A on Q' and Q 2 : 

KM = A(A), K-iM = A(A(e')), ^-a<A) = A(A(&")). 

We note atx> that oi = an , 6i = an and A r = {2}. This means that we can 
write (5.7.3) as 

A(A) = (A - a u ) A (A(^)) - £ < A (A (<?)). (5.7.4) 

><AT 

Now consider the star. The equation corresponding to (5.7.3) is equation 
(5.3.5): 

A-q l -g-iL = ffij^>. (5.7.5) 

IT7.7(a-p>) 

To rewrite this in the same notation as (5.7.4), we note that for a star on 
m + 1 = n vertices, with the centre labelled 1, M = {2,3,...,m + 1), ai = 
an, bi = 01,4+1, so that 

A(A) = <A-an)A(A{S’))-£ a ?r [] A(A(ei))- (5-7.6) 

Note that we have assigned the m = n — 1 p’s to the m connected components 
of Q’ so that p, is assigned to <71, i, i = 1,2 ,. ..,m. Note also that although 
the first tarns on the right of (5.7.4), (5.7.6) are identical, the second terms are 
different. For the star, the vertex set of Q" is empty. 

Par ter (1960) (265) obtained a general result which embraces the particular 
cases (5.7.4) and (5.7.6): 

Lemma 5.7.1 A(A) = (A-a 11 )A(A(«'))-5: >cAf a?,A(A(aj')).n« /A , A(A(<7;» 
itith the convention that A(A(<7")) = 1 if at for the star, the vertex set of Q, u empty. 

Lemma 5.7.1, like the corresponding result (5.7.5) for the star, is effectively 
a partial fraction expansion. In the general case, it is 
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where we haw used the fact that Q' has m separate connected components </J, 
so that 

A(A(G')> = I] A < A <3))- 

i<AT 

Equation (5.7.7) provides the basis for an inductive argument: deleting vertex 
1 of Q splits Q' into m components Q J, and bears the same relation to ^ as 
Q ' does to &. This means that if we can effect the reconstruction of A on the 
components Q 1 - of Q' from data referring to Q" and Q'-* then we can reconstruct 
the whole of A. 

Now* since 5' is itself a tree, and Q” is obtained by deleting vertex j from 
Gj f the roots of A(A((JJ')) should interlace the roots of A(A(£J)), just as the p- 
interlace the A*, i.e., (5.7.2). But equation (5.7.7) gives A(A(£p) as a result 
of the partial fraction expansion. We are given 

A{A) = /(A) = f[(A-A i ) ( (5.7.8) 

M 



A(A($')) = g( A) = jl(A - *). (5.7.9) 

M 

Now we must assign the n — 1 /i’s among the m components 5J. Suppose 
has Vj vertices then w r e must split the indices {1,2 ,.. ., n — 1} into m sets, so 
that if j € then Q J is assigned Vj indices. This is equivalent to grouping the 
terms in p(A) into m terms &(A), wrhere gj (A) has degree Vj: 

9( A) = n 



This means that we must check to see if, when /(A)/p(A) ts expanded into 
partial fractions, as 



/(A) . V' A,-(A) 



(5.7.10) 



where fy(A) is a monic polynomial with deg(fy) < deg(^), and if the A ? s and 
p’s interlace as in (5.7.2), then is positive, and the zeros of fy(A) and gj( A) 
interlace. To do this, it is best to change back into a form like Lemma 5.7.1 by 
multiplying throughout by <?(A): 



/<A) = (A- a ) S (A)-5:y J MA)« i (A) 



where ti,(A) = ^(A)/g,(A). We note that 



u i<A) = 

•GO 
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where Q consists of { 1, 2, . . . , n — 1} less those indices which are assigned to <fy. 
Choose j € A r and suppose that p r ,p r + p are two successive zeros of then, 
shrce g(p r ) = 0 = gUs rtft ), we have 

/(Mr) = -WMMr)Uj(Mr) 

/(Mr*p) = -yj/»j(Mr4 P )«7(/*r*p). 



We need to show that hj{ X) has a zero between p r and /J r<fpl t.e., that h J (/i r ) > 
(m„ p ) have opposite signs. The terms (/ i r - M.) and (M r4 „ - p, ) appearing in 
Uj (p r ) and ti^(/i r4|> ) will have the same sign except for those p t lying between 
p r4p and p r \ these are p — 1 such /i’s, with indices r + p — r + 2 ,r + 1 . 
Thus 

p odd ; Uj{u r ) y Uj(u v + p ) have the same sign 
p even ; Uj(u r ),Uj(u r + p ) have opposite signs. 

By assumption /(A) has just one zero between any two successive /i ? s; thus 



p odd ; /(p r ), /(/*„,,) have opposite signs 
p even ; /0i r ), /(/»„,,) have the same sign. 

Combining these results, we see that hj(p r ) and h ) lit rtp ) must have opposite 
signs. 

Now we check that y, is positive. Suppose Vj = q. and the roots of p>(A) 
and hj(A) are (p a ,)i and «*P®cti'ely, where 



M a , > > M ai > * • • > > /•«,,- 

Then 

, ,-i 

9jW = n <A ”*«.>• m a )= n^ A “*'-•) 

»-i t-i 

/(A) = f[(X-X<) 
i-1 



and suppose ai = r, so that 

/(M r ) = HIM, - A,). 

i-1 

Now* Ai, Aa, . . - , A r are all greater than p r , so that the sign of / (/v) is (-) r . All 
the v Qi are smaller than ^ so that hj^) > 0. Finally 

M/v)=n^-*) 

where the sum is taken over those i € {1,2, . .. ,n — l}/{a\ y ct2, . . . ,o„}. But 
for the sign we need to consider only those p, > p^.; there are r — 1 of these, so 
that the sign of Uj(p r ) is (— ) r “*. Thus y; > 0. 




116 



Chapter 5 



This yields the first stage in the construction of A: take /(A), 3(A) and form 
the partial fraction expansion (5.7.10); an = a, a 1; - = (/j) 1 ' 2 , while the zeros 
of 3,(A) and h,(A) are the eigenvalues of the components of A,j. 

Figure 5.7.1 shows an example of a tree. 



1 




Figure 5.7.4 - A tree on 9 vertices. 



The matrix has the form 



.4 = 



XX X 

X X X X X 
X X 
X X 
X X 

X XXX 

X X 

X XX 
X X 



In stage 1, A,i has two components; we find and we find new data 

which will allow us to construct the star on vertices {2, 3, 4, 5}, and the star- 
path on vertices {6, 7,8, 9}. To carry out the second stage we can, if we choose, 
relabel each of the connected components so that 2 — * 1 and 6 — * 1. 

We haw assumed that the data for constructing A is two strictly interlacing 
spectra. However, as with the path and the star, it is possible to use one 
spectrum <»(A) = (A,-)" and the first coefficients uu, i = 1,2, ...,n of the 
normalised eigenvectors of A. instead. We recall the result proved for a general 
matrix A 6 Sn, namely that the eigenvalues of A,i are the zeros of 






z- 

i-1 



A-A, 



= 0 . 



Further discussion of, and reference to, eigenvalue problems related to trees may 
be found in Nylen and Uhlig (1994) (252). 

Further references to the vast literature on inverse eigenvalue problems may 
be found in Gladwell (1986a) (107), GladweU (1996) (124), Nocedal and Over- 
ton (1983) (251), Friedland. Nocedal and Overton (1987) (95), Ikramov and 
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Chugunov (2000) |184|. Xu (1998) [339). Chu (1998) (58) and Chu and Golub 
(2001) (59). 



1. Show that if the eigenvalues A, of A € S n have maximum modulus A, then 
Kjl < A for alltj = 1,2,. ...n. 

2. Show that if A € S H is a matrix on Q. then the maximum number of (non- 
diagonal) zero entries in A is attained when Q is a tree, and is n 2 — 3n 4- 2. 

3. Construct an algorithm to form A from (A*)”, given the structure 
of Q. Use it to construct A on the graph Q of Figure 5.7.4. Take 
{A*}? = <1,3, 5, 7, 9, 11, 13, 15, 17}, {/r,}f = {2, 4,6, 8, 10, 12,14, 16}. As 
a check, find the eigenvalues of A and Ai. 




Chapter 6 

Positivity 



There are then two kinds of intellect: the one able to penetrate acutely and 
deeply into the conclusions of given premises, and this is the precise intellect; 
the other able to comprehend a great number of premises without confusing 

them, and this is the mathematical intellect. 

Pascal’s Penstes, 2 



6.1 Introduction 

The basic eigenvalue analysis of real symmetric matrices was discussed in Chap- 
ter 1. The eigenvalue properties described there are shared by all positive- 
definite (or semi-definite) matrices. This Chapter, which may be missed on a 
first reading, provides proofs of some of the results which were used in Chap- 
ter 1. Foremost among these are Theorem 6.3.1, that if A € S*, then it has 
n real eigenvectors which are orthonormal, and thus span V u \ and Theorem 
6.3.7 that provides necessary' and sufficient conditions for the matrix A to be 
positive-definite. Signs, positive or negative, provide the recurring theme for 
this Chapter, and hence our choice for the Chapter heading: posithity. 

In Chapter 3 we focussed our attention on a narrower class, Jacobi ma- 
trices, and found that they had additional eigen-properties: they had distinct 
eigenvalues and, with increasing t, the eigenvector u, became increasingly os- 
cillatory', meaning that there w r as an increasing number of sign changes among 
the elements It will be showm in this Chapter that many 

of the eigen-properties of such matrices are shared by a wider class of so- 
called oscillatory matrices. Actually, there are twin classes of matrices, os- 
cillatory' and sign- oscillatory, as described in Section 6.5. If A is oscillatory', 
and Z = diag{ 1,-1, l,...,(-) n " 1 ), then A = ZAZ is sign-oscillatory', and vice 
versa. The Jacobi matrix J of equation (3.1.4) is actually sign-oscillatory. 
These matrices were introduced and extensively studied by Gantmacher and 
Krein (1950) [98], see also Gantmacher (1959) (97). The matrices appearing in 
lumped-maas or finite element models of strings, rods and beams are all oscilla- 
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tory or sign-oscillatory; this Chapter serves as reference material for the study 
of oscillatory matrices. 

The theorem upon which the whole of the analysis of oscillatory matrices 
depends, is Perron’s theorem (Theorem 6.5.1). This relates to a strictly positive 
matrix, one that has all its elements strictly positive, and states that such a 
matrix has one eigenvalue, the greatest in magnitude, that is real and positive; 
the corresponding eigenvector has all its coefficients strictly positive. 

The matrices appearing in mechanics are usually not strictly positive; such 
matrices appear in Economics and Operational Analysis. Instead, the matrices 
are oscillatory. (See the precise definition in Section 6.6.1.) In order to apply 
Perron’s theorem to such matrices, we need two essential steps. First, if A 
is oscillatory, then B = A n_I is totally positive (TP). This term, which is 
introduced in Section 6.6.1, means that not only all the elements of B are strictly 
positive, but also all the minor* (Section 6.2) of B. Note that the eigenvalues of 
B are the (n — l)th powers, , of the eigenvalues of A, while its eigenvectors 
are the eigenvectors of A. The other step that is needed is the introduction 
of the concept of a compound matrix (Section 6.2). The compound matrix 
A„ is formed from all the N = (JJ) pth-order minors of A. The important 
Binet-Cauchy Theorem, Theorem 6.2.3, shows (Ex. 6.3.1) that the eigenvalues 
of Ay are simply products of p eigenvalues of A. The argument then runs 
as follows. Suppose A is oscillatory, then B = A" -1 is TP, and hence, for 
p = 1,2, . . . ,n, B P is strictly positive (not TP). The first conclusion (Theorem 
6.10.1) is that the eigenvalues of A are positive and distinct, like those of J or 

Before beginning the analysis proper, we point out a notational matter which 
must be understood if confusion is to be avoided. In Chapter 3, in dealing 
with a Jacobi matrix J, a positive semi-definite tridiagonal matrix with negative 
codiagonal, the eigenvalues were labelled in increasing order, i.e., 0 < Ai < A 2 < 
••• < A*. The eigenvectors then became increasingly oscillatory, as described 
in Theorem 3.3.1. In Ex. 3.3.2, it was pointed out that if the eigenvalues of 
J = ZJZ. a positive semi-definite tridiagonal matrix with positive codiagonal 
(an oscillatory matrix if it is actually non-singular, i.e., positive-definite) are 
labelled in decreasing order, i.e., Ai > A 2 > ■ - ■ > A„ > 0, then the eigenvectors 
still satisfy Theorem 3.3.1. In this Chapter, in dealing with oscillatory matrices, 
we shall keep the same ordering, i.e., A, > A a > • - • > A„ > 0. Theorem 6.10.2 
is a generalisation of Theorem 3.3.1. 



6.2 Minors 



Suppose A € A/ n . To gain some insight into the structure of A, and into the 
relative sizes of its elements, we introduce the concept of a minor. A minor of 
order p of the matrix A is the determinant constructed from the elements of A 
in p different rows and p different columns. Thus, the elements of A themselves 
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are minors of order 1, while det(A) is the only minor of order n; 013 , 
and det( A) are all minors of A. 

Folkwing Ando (1987) [4] we let Q p%n% with 1 < p < n, denote the set 
of strictly increasing sequences a of p integers a lt a 2 ? ... t a^ taken horn cj = 
{l,2,...,n}. The complement c/of a is the increasingly arranged sequence 
{1,2,..., n}\a = u/\a, so that o' € Q„- V . n . When a € Q p .„, 0 € Q q , n and 
a n 0 = 0, their union, a U p, should always be rearranged increasingly to 
become an element of Q Ttn {r = P*M). We will often use two special sequences: 
0 = 0(p) = (1,2, . - . ,p} and = <Xp) = {n-p+ 1, . . . , n}, and their complements 
^ = ^(p) = (p + 1, . . . , n}, = tf>'(p) = ( l t 2 , . . . , n - p}. When the argument 

is omitted in 0 or & it will be understood to be p. 

The submatrix formed from rows a and columns p of A is denoted by A[a|0); 
A[a|o) is written A[a]. The minor of A taken from rows a and columns p is 
denoted by A{a\0)\ thus 






“a, A 


a -A, • 




A{a,0) = 


“a,*, 


a “jAi - 


■ a oaA» 




SA 




■ 



( 6 . 2 . 1 ) 



The minor -4(a;<*) is abbreviated to A(a). 

The cofactor of a^, introduced in Section 1.3, is a minor with a sign attached 
to it: 

A* = ( 6 - 2 . 2 ) 

where 

(6.2.3) 

andt' = {l, 2 ,...,«-l,» + l,...,n}=w\i, ? = {l, 2 ,...,i- l,j + l,...,n} = 

<Jj\ &ij is sometimes called the minor of a^. 

If A € A/*,, then we can form a new r matrix A = (d*/) from the minors of 
elements of A. We may prow 



Theorem 6.2.1 Let A = (d tJ ), then the minors of A are given by 



i(a;0) = (det(A)r'.4(a';/r). 



Proof. Consider the theorem for a = 0 = 0; the general case may be 
obtained by a suitable rearrangement of the rows and columns. Since dij = 
(-) 4+) 'A> . we may write 



A n Ai2 
A 2 \ A 22 

A V < 



, ,p 

A *p 

■V 



( 6 . 2 , 1 ) 



B = A(a,0) = 
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Multiplying this by det(A), and writing the determinant in (6.2.4) as that of an 
n x n matrix, we find 





f 


’ *11 
*ai 


*12 . 
*22 - 


A\p Ai.r+x • 
. *2.p+i 


. *ln 
• *2« 


B ■ det(A) =det 




Y 


V: 


: v r > : 


• 

. 0 




v 


. • 


o .. 


. 0 0 


. 1 



*11 


*13 • 


• *!• " 






*12 * 


• *2., 




. *nl 


- 


• *•« - 













so that on using equation (1.3.10) we obtain 



B • det(A) 



det(A) 


0 


0 


0 


0 


det(A) ... 


0 


0 




det(A) 


0 


0 




a2j.*i 


a^n^.i 


... On.pt I 


ain 


02,, 




- Inn 



(det(A 



so that the theorem holds when det(A) ^ 0. Continuity considerations show 
that the theorem also holds when det(A) = 0. B 

One of the implications of this theorem is that when det(A) = 0 t the rank 
of A is at most l t meaning that all the row’s of A are multiples of each other, 
as are all the columns. There is another corollary 



Corollary 6.2.1 

det{A) = (det(A))"- 1 . 



There is another way to form a matrix from minors of a given matrix. Sup- 
pose A € M„ and 1< p < n, and put a = S(p) := {1, 2, . . . ,p}. We can define 
hi, by 

by = A(0U«, SU j), i,j =p+l,p + 2,...,n. 



The matrix B 6 M n -p. Thus, if p = 2 and 




2 

1 

1 

0 



3 
-1 

4 
3 



4 

2 

1 

2 



-5 

-2 




The matrix B is called a bordered matrix; the indices i, j border 0. 
Sylvester’s Theorem on bordered determinants is 



Theorem 6.2.2 Suppose A € Af„, 1 < P < n, fey = A{9 U t, 0 U j) /or 
i,j = p + l,p 4- 2 n, and B = (by), then 

det(B) = = (A(f>;fl))"-^‘d«(A). 
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Proof. Theorem 6.2.1, with p replaced by n -p- 1, and a = {p+l,...,r- 
l,r + l n) = 0 '\r, 0 = {p + 1 ,. . .,a - l,s + 1 ,.. .,n} = *V, shows that 



•V. = i( a ;^) = (det(A))"-'- 2 A(flUr;eus) 
= (det(A))-»- a h r .. 


(6.2.5) 


The Corollary of Theorem 6.2.1 shows that if C = (c,.)^,,, then 




det(C) = (i^y))-^*. 


(6.2.6) 


But according to (6.2.5), 




det(C) = B(^;0')(det(A)) < "- , ’- 2,( "-' , 
= det(B)(det(A)) ( "-'- aMn -'’> 


(6.2.7) 


and from Theorem 6.2.1 




A{O':0') = (det(A)) n -'- 1 A{e;O) 


(6.2.8) 



so that on substituting (6.2.8) into (6.2.6) we find 

det(C) = (det(A)) < "“' , - 1>> (^(^;6»))"-»’- 1 

which, on comparison with (6.2.7), yields the required result when det(A) ^ 0. 
(Continuity considerations show' that the theorem still holds when det(A) = 0. 

■ 

Corollary 6.2.2 // a,0 € 0fla = O. 0n0 = O, then 
Bfr0) = (A{9;O)r-'A(O Uo;0U0). 

Corollary 6.2.3 Suppose € Qpjt ani * 

b ti = A(aUi,0Uj) 

fori =<tp + j = 0 p + l,...,n and y,S € Q 9M with -y, > a,, 6, > 0 p , 

then 

B(r,S) = (A(a;ft))*- l A(aUr,0U6)- 

This is the general form of Sylvester’s Theorem. For a proof, see 
Gantmadxer (1959) (97], Vol. I, p. 32. 

We now introduce the powerful Binet- Cauchy Theorem. 

Theorem 6.2.3 // A 6 M mJ) , B € A/*.„ and C = AB. a € Q,, mt 0 € Q pn 
then 

C(o,0) = Y J M«^)B{T,0) (6.2.9) 

utftere the sum extends to all 7 € Qp,k- 
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The theorem is a generalisation of the formula for namely 

i 

^ = EVw- 

p-i 

The proof may be found in Gantmacher (1959) [971 . VoL I, p. 9. 

The importance of the Binet-Cauchy Theorem lies in its application to com- 
pound matrices, which we now define. 

Suppose first that A is square, i.e., A € A/*. We shall define the compound 
matrix A v . Consider all the sequences o € Qp % n\ there are 

w ■ G) " 

such sequences. For given n,p, ihe A’ sequences may be arranged in ascending 
order 1,2, . ..,JV. This may be done by associating with the sequence a = 
{o 1 ,i» 2 ,. .. ,a t ) the number with digits a 1( a 2 , ... ,o v in the base of d = N + 1 . 
This procedure associates a specific index s = s(o) with each sequence o; s lies 
in the range 1 < s < N. Thus, when n = 5, p = 3, we have iV = 10, and 
the combinations are 123, 124, 125, 134, 135, 145, 234, 235, 245, 315. Thus 
*(124) = 2, while s(245) = 9. The element a,, of Ay is then given by 

a,, = A(a y 0) 

where s = s(a), I = s(0). 

Although we shall not need it in this book, a compound matrix can be defined 
for a rectangular matrix A € Now 

6 JfuN ' w = (p)’ A ' = C) 

and a,, is given by (6.2.9) for a € Q v , m . P 6 Q p ,„. The Binet-Cauchy Theorem 
now” states 



Theorem 6.2.4 // A 6 B € C = AB andp < min (m,A;,n), then 

Cp = A^e- 

Proof. The equation (6.2.9) may be written 

* 

•V. = X] a ^ 6 '* 

r-i 

where r = s(a), * = s(p), t = s(y). m 

Corollary 6.2.4 If A € Af„ is non- singular then the pth order compound ma- 
trix of A -1 is the m terse of the pth order compound matrix of A. 

Proof. Let B = A -1 , then AB = I implies A,.B V = I? so that B ; = 

(A)-'- ■ 




124 



Chapter 6 



Exercises 6.2 

1. If A 6 M,i is non-singular, then equation (1.3.20) shows that its inverse 
R = A -1 has elements 

r *i = Ajif det(A) = (-) <+ ^,/det(A). 

Use Theorem 6.2.1 to show that if a, 3 € Q p .„ then 

det(A)R(a;fl) = (-)'A(o';^) 

where 

f= £(o m +fl m ). 

m— 1 

2. 110= {1,2,... ,p) and 0 = {n -p+ 1,.. ,,n} then A{0,<j>) and A(frB) are 
called the pth order comer minors of A. Use Ex. 6.2.1 to show that the 
corner minors of R are given by 

det(A) • fi(M) = 

where i = r(p -f 1). Note that A{&\ <*/) is an (n — p)th order corner minor 
of A. 

3. Equations (1.3.7), (1.3.8) are a particular case of laplace's expansion of a 
determinant t 

det{A) = 

where a € Qp % u is fixed, the sum is taken over all 3 € Qp % n and t = 
(a m + 0 W ). Establish this result and show that there is a similar 
expansion with fixed and a varying over Q v .n. 

4. Suppose A € M r . Use the Binet-Cauchy theorem to show that the pth 

compound matrix of A m is i.e., 

(A*) P = (A P r=A?. 

5. Use the Binet-Cauchy theorem to show* that if Q is an orthogonal matrix, 
then so is Gp, the pth compound matrix of Q. 

6. If B = A" 4 , write the minors of B in terms of the minors of A; use the 
notation (6.2.9) to show that 

B( Q J) = Y, y)Mr, V) - • • -4(y ,m -"i P) 

where the sum is over all 7; . . . y <m-1) € Q p .n- 
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6.3 A general representation of a symmetric ma- 
trix 

We begin with two theorems. 

Theorem 6.3.1 If Af S n , then A has n real eigenvectors forming on ortho- 
normal system. 

Theorem 6.3.2 To each m-fold eigenvalue A 0 of A € S n , there correspond tn 
linearly independent eigenvectors. 

In Section 1.4 we showed that the eigenvectors corresponding to distinct 
eigenvalues are orthogonal. This means that if all the eigenvalues of A are 
distinct, then it has n orthogonal eigenvectors which may be scaled so that they 
are orthonormal. It sufficies to prove Theorem 6.3.2. 

Proof. Suppose that A 0 is an m-fold eigenvalue of A, i.e., A(A) = det(A-AI) 
has Ao as an m-fold root, and that B = A — AqI has rank p, so that the equation 

Bu = (A — Aol)u = 0 (6.3.1) 

has r = n —p linearly independent solutions. We need to prove that r = m. 
Now 

A(A) = det(A - AI) = det(B - (A - Ao)I) 

= £?^<-)%(A-Ao)<, i = n - t, 

where Tj is the sum of the jth-order principal minors of B, and To = 1. But B 
has rank p, so that T„ = 0 = T*_i = •*• = T v + \ and therefore 



ao that m > r. It ts sufficient to prow that T p ^ 0, for then A(A) will have a 
r-fold root, i.e. t m = r. Without loss of generality we may assume that the first 
p rows of B are linearly independent, so that any row of B may be expressed as 
a linear combination of the first p rows, i.e., 

p 

bij = *.j = 1,2, ...,n, 

*-l 



which may be written 



B = CB 0 . 



(6.3.2) 



where C € M UJ ,, and Bo 6 A/p,* is formed from the first p rows of B. Now 
apply the Binet-Cauchy Theorem 6.2.3 to (6.3.2): 



B(a;0) = £C(o;q)Bo(->;0). (6.3.3) 



But C has only p columns, and similarly Bo has only p rows, and they are the 
rows 1,2, ... .p of B. Thus, there is only one term in the expansion (6.3.3): 



B(a;0) = C(a;*)B(0;£), (6.3.4) 
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where 0 = {1, 2, . . . ,p}. Similarly 

Btf;0) = C(0,O)B(O;0}. (6.3.5) 

But B b symmetric, so that B(3;0) = B(0:0), and thus, on combining (6.3.4), 
(6.3.5), we haw 

B(a; 9) = C(a\O)C{0\O)Bi,0,O). (6.3.6) 

All the minors on the left cannot vanish, since then B would haw rank leas than 
p; we must have B(0;0) ^ 0. But then (6.3.6) with 3 = c> gives 

B(a; a )=(C(«*;0)) a B(M). 

This means that all the pi h order principal minors of B have the same sign, and 
one at least, B(0;0) is non-zero. Thus T Vy their sum, must be non-zero. Hence 
m = r. ■ 

We may now assert that if A € S*, then it has n eigenvalues (A ,)" and n 
orthononnal eigenvectors (u»)j . This means that 

Au» = A^, t = l,2,...,n, 

which may be combined to yield 

AU = U A, UU T = U r U = I (6.3.7) 

and this may be transformed to 

A = UAU t . (6.3.8) 

This is a most important representation of a symmetric matrix. 



Exercises 6.3 



1. Apply the Binet-Cauchy Theorem, in the form of Theorem 6.2.4, to equa- 
tion (6.3.7) to show that the eigenvalues of A? are all the products A^A,^ 

-v 

2. Show that the eigenvectors of Ay are the columns of the compound matrix 

u p . 



6.4 Quadratic forms 

Suppose A € Sn, then 

A(x,x) = x r Ax = a„xj + 2a I2 x,x 2 + • •• + 2a„. n _ 1 x,,_ 1 x n + (6.4.1) 

is called the quadratic farm associated with A. One of our aims in this section 
is to find necessary and sufficient conditions for A to be positive definite (PD) f 
i.e., A(x.x) > 0 for all x ^ 0. 
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First, we consider a number of different ways of expressing A{x, x]. Let 



n 

M x ) = 52 a it x t> * = 1,2, ... t n, 

i - > 




This yields 



A(xx) = £a^(x). 

»-> 



Oil 


o 12 


Olw 


Ai(*) 


021 


022 


02» 


A»(*) 


o»i 


On 2 


• Onn 


A n {x) 


Ailx) 


A 2 (x) 


■ A n (x) 


M*,*) 



since the last column is a combination of the first n columns, and 



( 6 , 1 . 2 ) 



(6,1.3) 



( 6 , 1 , 1 ) 



Theorem 6.4.1 // det(A) ^ 0, then 





Oil 


012 


Oln 


Mx) 


-1 


021 


d22 ... 


02 n 


A 2 {x) 


det(A) 


Onl 


On 2 


Onn 


A„(x) 




Ai(x) 


A 2 {x) ... 


A n (x) 


0 



(6,1.5) 



Proof. Expand the zero determinant (6.4.4) along its last row. ■ 
Now* we introduce the quantities 



X,(x) = ^ 1 (x),X 2 (x) = 



“11 

“21 



etc., up to AT r( (x), and prove 





“11 


“12 


A,(X) 


,* 8 (x) = 


“21 


“22 


A 2 (x) 




“31 


“32 


A*(x) 



( 6 , 1 . 6 ) 



Theorem 6.4.2 // 0 = {l,2,...,p} and D v = A(0:O) ? 0, p = 1,2,..., n 
then the (X,,(x))J are linearly Independent. 

Proof. We see that Xi(x) = X^-i a 'i x i< while X 2 (x) = £"_ 2 (ana2> - 
ana\j)xj, and generally 

*p(x) = D v x p + terms in j fM1 ,...,x n . 

Thus we see that in the reversed sequence A'„(x), A' n -i(x), . . . , Xi(x), each term 
involves one more Xj than the previous one. This means that the (AT*(x))J can 
all be simultaneously zero iff all the (x,)“ are 2 ero. ■ 

This leads us to an important expression for A(x, x) given by 
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Theorem 6.4.3 (Jacobi). 
0, p = 1,2 ,...,n then 



If D 0 = 1, 0 = {1,2,..., p} and D„ = A{9,9) ji 






r-i 



<*,(*))* 

DpDp \ ' 



(6,1.7) 



Note that, on account of Theorem 6.4.2, this equation expresses A(x, x) as 
a sum of multiples of squares of linearly independent combinations of the (a*)". 
Proof. Put ft = 0. and 



<*n 


a l2 -• 


"I* 


A,(x) 


«J1 


<*22 -• 




A 2 (x) 


a pl 


<*v2 -• 


Opp 


A y (x) 


Mix) 


A 2 (x) .. 


Ap(x) 


0 



( 6 , 1 . 8 ) 



and find the recurrence relation linking P p and P p _i . P p (x, x) is the determinant 
of a symmetric matrix C € Sp+i, i.e., 



P »(x. x) = C(0(p 4- 1); 0[p 4- 1)). 



Apply Theorem 6.2.2 to this, letting 



fty = C(9{p - 1) U i; 9(p - 1) Uj), «, j = p,p + 1 



del(B) = byyby.l'fH — Vpt l^Pt l.»' , fl Q, 

= C(9(P-1);9(P-1)).C(«( P + 1);«(P+1)) 

But b„ = D„, 6 pl i. p „ = P p _i(x,x), 6 p . p _, = 6 p _,. p = AT p (x) whUe C(0(p- 
l)i 0(p - 1)) = D„-i, C(9(p+ l);0(p + 1)) = /^Ix.x). Thus, equation (6,1.9) 

D p P p -.(x,x) - Xl(x) = D p .,P p (x,x) 
or, since the D v are non-zero 



-p p M = -/y-,(x,x) + p= 1[2> ... tn . 



D. 



O p .i D p D p -i 



( 6 , 1 . 10 ) 



Now Theorem (6,1.1) states that 



A(x,x) = 



-P,,(x,x) 

A." 



so that on summing equation (6,1.10) from 1 to n and using P 0 = 0 we find the 
required sum (6.4.7). ■ 



Theorem 6.4.4 Suppose A € S*, then A is PD iff D, > 0, t = 1,2, . .. ,n. 
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Proof. First we prove that if A is PD, then det(A) > 0. Since A € S n , it 
has, by the Corollary to Theorem 6.3.2, n eigenvalues (Aj)f and n oithonormal 
eigenvectors (u,-)y such that Au, = A.u,. Thus A* = ufAu, = A(u,,u,) > 0 
and therefore det(A) = JJ"., A,- > 0, i.e., D n > 0. 

If A is PD, then the matrix obtained by deleting the last ) rows and columns 
of A is PD, for j = 1, 1 ,. . . ,n — 1. Therefore, their determinants are positive, 
i.e., > 0. We have proved that if A is PD, then (D,)" > 0. 

Now suppose that (D,)J > 0, then equation (6.4.7) shows that A(x,x) > 0, 
for the (Af,(x))i can be simultaneously zero only if x = 0. Thus A is PD. ■ 

Corollary 6.4.1 If A € Sn is PD, then alt the principal minors A(o;o) = 
A(o), a € Qp.,„ p = 1,2, . . . ,n, are positive. 

If A € S H is merely positive semi-definite (PSD), then the leading principal 
minors, and indeed all the principal minors are non-negative. We prove 

Theorem 6.4.5 If A € Sn it PSD and, for some p satisfying 1 < p < n, D p = 
A(0; 0) = 0, then every principal minor bordering D v is zero. In particular, the 
leading principal minors R,, p < q < n, are zero. 

We prove that the D q are zero, and leave the remaining result to an Exercise. 
Proof. There are two cases: 

i) p = 1, then Di = an = 0, and 

<*i> 

a H a a 

implies (ay)" = 0, so that (D v )" = 0; in this case xi does not appear in .4(x,x) 
at all. 

ii) an 0 and, for some p, 1 < p < n — 1, D v 0, D p * i =0. (If p = n — 1, 
there is nothing further to prove; we may therefore take p < n — 1.) 

We introduce bordered determinants 

b tj = A(0U.;0Uj), «,j=p+l n 

and form B = (fey)",,. By Sylvester’s identity (Corollary to Theorem 6.2.2), 
if Qi > p and o € Q,, n < r < n — p, then 

B( o;o) = (4(#;8)rM(6Ua;flUa) 

so that B is PSD. Since bp.ip.i = D pJ1 = 0, the matrix falls under case 1 and 
if q >p + 1, r = q-p- 1, a = {p+ 1,...,$}, then 

0 = B(a;a) = {A(0{p)Ap))Y A{O{q)-0{q)) 

so that A(0{p)\ fl(p)) = D p 0 implies A(0(g);0(g)) = D q = 0. ■ 

This theorem implies that if A € Sn is PSD. then, for ®me p. 1 < p < 

n, (R)i > 0, (D,)J 41 =0. 
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Exercises 6.4 

1. Show that A € S„ is PD iff its eigenvalues (A,)J are positive; it is PSD iff 
its eigenvalues are non-negative. 

2. Show that if A e S„ is PSD then A is singular, and that x r Ax = 0 iff 
Ax = 0. 

3. Show that if A e S„ is PSD and if a ti = 0 for some « satisfying 1 < i < n, 
then a,, = 0 for j = 1, 2, . . . , n. This means that if a,, = 0 then i, does 
not appear in x r Ax. 

4. Show that if A 6 S n is PSD and has rank r then it has a positive principal 
minor of order r. 

These examples are merely a selection of properties of PD and PSD matrices 
to be found in Chapter 7 of Horn and Johnson (1985) (183). 

6.5 Perron’s theorem 

Mast matrices appearing in classical vibration problems are symmetric. It is 
therefore known that they have real eigenvalues, and a complete set of ortho- 
normal eigenvectors. Often the matrices are PD, so that their eigenvalues, in 
addition to being real, are positive. However, the whole theory relating to os- 
cillatory matrices depends on a basic result relating to a class of not necessarily 
symmetric matrices, as we now describe. 

We recall a>me definitions. If a vector x has all its elements positive (non- 
negative) we shall say x > 0 (> 0) and shall say that x is positive (non -negative). 
If x, y are in V* then x > y is equivalent to x — y > 0. The matrix A £ M r4 is 
said to be positive (non-negariue) if <Hj > 0 (> 0) for ail i,j = 1, 2,. . . ,n. 

Up to this point the only norm we have used for a vector x € Vn is the 
Euclidean, or so-called I? norm: 

11*112 = <][>| 2 )'. <6.5.1) 

i-1 

We can define the L 2 norm of a matrix A € M„: 

l|A||a = <EMV. (6.5.2) 

The norm is variously called the Frobenius norm. Schur norm or Hilbert-Schmxdt 
norm. 

We will need another norm, the L\ norm: 



INI* = £m, 



<6.5.3) 
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II A||» = £ |oy|. (6.5.4) 

A norm is like a distance, as such it must satisfy various fundamental conditions, 
for which see Ex. 6.5.1. For a definitive and extensive study of vector and matrix 
norms, see Horn and Johnson (1985) (183], Section 5.6. 

We may now prove Perron's theorem, following Bellman (1970) (25]. 

Theorem 6.5.1 (Perron). Suppose A € M„ and A > 0. Then A has 
a unique eigenvalue p which has greatest absolute value. This eigenvalue is 
positive and simple , and its associated eigenvector can be taken to be positive. 
The eigenvalue p is often called the Perron root of A. 

Proof. Let 5(A) be the set of all non-negative A for which there exist non- 
negative x such that 

Ax > Ax (6.5.5) 

We shall consider only Li-normalised vectors x, i.e., such that ||x||i = z, = 
1. (Since x > 0, |zi| = z,.) This therefore excludes the zero vector. If x 
satisfies (6.5.5), then Ex. 6.5.2 shows that 

A||x||, < ||Ax||, < || A||i.||x||i, (6.5.6) 

so that 0 < A < HAH,. This shows that the set 5(A) is bounded. It is clearly 
not empty, because A is positive. The bounded set 5(A) has a least upper 
bound; let it be A<>. Let A 1( A 2 , ... be a sequence of A’s in 5(A) converging to 
Ao, and x, a corresponding sequence of x’s satisfying Ax, > A,x,-. The set of 
all x such that ||x||) = 1 is closed and bounded; therefore, the sequence (x,) 
contains a convergent sequence {x„ ( } converging to a non-negative vector Xq 
satisfying ||xo||, = 1, and (Ex. 6.5.3) 

Axo > Aoxo. (6.5.7) 

This means that Ao € 5(A). We shall now show that equality holds in (6.5.7), 
and we do so by reduction to a contradiction. 

Let 

d = Axo — AoXo > 0 
and suppose one of the d,. say dj, is positive. Put 

Vi = r*t + (d*/2AoMii 

then the «th row of Ay — Aoy is 

e,- = di + a„d J /(2A 0 ) - dfi^f 2 > 0. 

Now let A = Ao + min, •(«>,-/&), then A > Ao, and 

Ay -Ay = e-(A-Ao)y 

= e - min,(e,/p,)y > 0 ' 
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This slates that A € 5(A), and that A is greater than the least upper bound, Aq, 
of 5(A). This contradiction implies that there is equality in equation (6.5.7), 

Axo = AoXo. (6.5.8) 

Thus is an eigenvalue and x 0 is an eigenvector, and Xq is necessarily positive 
(Ex. 6.5.4). We will show that Aq is the required Perron root. 

Suppose that there is another eigenvalue A, possibly complex, such that |A| > 
Ao, with z 0 being an associated eigenvector, so that Az = Az. Let |z| denote 
the vector with elements |zi|, |za|,. . • , |r.,|, then we deduce that 

|A| |z| = |Az| < A|z|. (6.5.9) 



But then the maximum property of Ao implies |A| < Ao, and hence |A| = Ao. 
Now the argument used earlier shows that equality holds in equation (6.5.9), 



But then 



A|z| = Ao|z|, 
|Az| = A|: 



z|>0. 



(6.5.10) 



and (Ex. 6.5.5) this can hold only if z = <w, where c is complex and w is 
positive; and this implies that A is positive, i.e., A = A^. We now show that Xq 
and w, both positive and both eigenvectors corresponding to A<>, are equivalent. 
Put y = x, — iw. and take 



e = uun(r,o/V',) = Xp/wj, 

then y is a non-negative eigenvector corresponding to Aq with y, = 0, so that 

+ *j7Vl + • • • + <*jnVn = 0, 

and since a$i > 0 for i = 1, 2, ... ,n we must have y = 0. Thus Xo = tw so 
that Ao is a simple eigenvalue. Thus Ao has all the properties asserted for the 
Perron root p. ■ 



Exercises 6.5 

1. A vector norm must satisfy three conditions: 

a) x|| > 0, and ||x|| = 0 iff x = 0 

b) bc|| = M • INI 

c) x + y||< ||*|| + IMI 

Show that both the L\ and the L 2 norm satisfy these conditions. 

2. Show that A € A/*, x € then 



IIAx|| 1 <||A|| I .||x|| 1 . 
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3. Verify that the vector Xo will in fact satisfy the inequality (6.5.7). 

4. Show that if x is a non-negative eigenvector of a positive matrix A € M„, 
then x > 0. This has the following logical negative consequence: 

if A > 0, x > 0, Ax = Ax 
and x, = 0 for some t = 1, . . . ,n, then x = 0. 

5. Show that if A € Af n is positive, then equation (6.5.10) can hold only if 
z = <w, where c is complex and w > 0. 



6.6 Totally non-negative matrices 

Suppose A € A/ n . The matrix A is said to be positive (see Section 6.5), written 
A > 0, if a,, > 0 for all i,j = l,2,...,n. Total positivity concerns all the 
minors of A, (see Section 6.2) not just its elements. If A € Af m . n , we say that 
A is 

1. TJV (totally non-negative) if alt the minors of A are non-negative; 

If A € M n , we say that A is 

2. NTN (non- singular and totally non-negative ) if A is non-singular and 

TN ; 

3. TP (totally positive) if all the minors are (strictly) positive: 

4. O (oscillatory) if A is TN, and a pouter, A m , is TP. 

Note that some authors, including ourselves in Glad well (1986b) |108), use 
totally positive (TP) instead of totally non-negative (TN), and strictly totally 
positive (STP) instead of totally positive (TP). Also, in Glad well (1986b) (108), 
following Gantmacher and Krein (1950) (98), we used completely instead of to- 
tally ; completely positive now has a quite different connotation. Reader, beware 
of these subtle distinctions! 

The concept of an oscillatory (or oscillation) matrix was effectively introduced 
by Gantmakher and Krein in the 1930’s, see Gantmacher and Krein (1950) (98). 
It was developed further by Gantmacher (1959) (97). The concept of total 
positivity had arisen much earlier than this, e.g., Fekete (1913) (86); it was 
first systematically explored by Karlin (1968) (190) in his book Total Positivity, 
Volume 1. (Volume 2 has never appeared!) Ando (1987) (4) reviews its history 
and proves important new results. All the concepts, total positivity, oscillatory, 
etc., arise in the study of in-line systems, rods, beams, splines, Sturm-Liouville 
differential equations, etc. 

The study of total positivity involves the delicate treatment of inequalities. 
Here are two typical examples, which the reader may verify: 
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i) if a > 0 or d > 0; 5 2 0 and c > 0; and 



a b 
c d 



> 0 , 



then a > 0 and d > 0; 



ii) if a > 0 or d > 0; b > 0 and c > 0; and 



a 

c 




then a > 0 and d > 0. 



The concept of total positivity is similar to positive-definiteness, but there 
are important differences between the two concepts: positive definiteness applies 
only to symmetric matrices. TP applies to any matrices in A/„; the condition 
for positive-definiteness involves only the principal minors, while TP involves all 
the minors. Clearly, if A 6 S„ is TN then it is PSD; if it is TP then it is PD; 
but the converses of these results are false. (Ex. 6.6.1). There is a theorem 
like Theorem 6.4.5 for matrices which are TN: 



Theorem 6.6.1 //A € M n , and A is TN , and A has a zero principal minor, 
then every minor bordering it is also zero. 



Proof. For simplicity we confine our attention to the leading principal mi- 
nors; this restriction can be removed at the expense of some complication in the 
argument. As in Theorem 6.1.5, there are two cases: 

1) D, = a u = 0. We assert that this implies that either (a,-,)y = 0 or 
( a ij)? =0- this were not true, then we could find a,! > 0 and a X j > 0 for 
&>me i,j satisfying 2^i <n, 2 < j < n. But then 



Oil 



a H 

a ’i 



= < 0 , 



which contradicts the statement that A is TN. Thus if a„ = 0 then either the 
first row of A or the first column of A must be zero. (See also Ex. 6.6.2.) 

2) a u ^ 0. Then for some p(l < p < n - 1) we have D v ^ 0, D v , x = 0. 
(Again, if p = n — 1, there is nothing further to prove.) We introduce bordered 
determinants 

bij = A(0Ui;0Uj) t,j' = p+ 



and form the matrix B = (6,*)". t . By Sylvester’s identity' (Corollary 6.2.3), if 
o,0 € Qr.n, Qi > P, 0i > P, then 



I}{a,tl) = {A(e;O)r- l A{9Va;0U0) 



so that B is TN. Since i = Dp.i = 0, the matrix falls under case 1. If 

9 > P + 1. and a = 0 = {p+ 1,.. then 

B(o;a) = <A(fl(p);6>(p)))*-»- , A(fl(9);fl(9)) = 0. 



But since D p = A{0{p),9{p)) 0, we have D q = 0. ■ 
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Theorem 6.6.2 If A 6 iV/„ and A is NTN, then all its principal minor* are 
positive. 



Proof. If any principal minor were zero, then, by Theorem 6.6.1, D„ = 
det(A) would be zero, but A is non-singular, so that det(A) ^ 0 (in fact 
det(A) > 0). ■ 

Corollary 6.6.1 If A € S n is NTH, then A is PD. 



If A € Af n is NTN, then we know that some elements of A are strictly 
positive, in particular, by Theorem 6.6.2, the diagonal elements <*„. We now 
prove an important result which shows that A will haw a so-called staircase 
structure. We first introduce some definitions. 

Let p= (pi,Pa,.-.,pn} be a sequence of integers bom {l,2,...,t»}. Then 
p is a staircase sequence if Pi < P2 < • • • < p„ and p t > i for all i = 1,2, ... ,n. 
Thus p = {2,3, 3,5, 5} is a staircase sequence. 

Suppose p and q are staircase sequences. A matrix A € M„ is said to 
be a p,q- staircase matrix if a tj = 0 when j > p, or i > qj. Suppose p = 
{2, 3, 3, 5, 5}, q = {2,4,5, 5, 5} then 



**n 

<*2i 



a, 2 
<*22 

<*32 



<*42 



<*23 

<*43 <*44 

<*63 <*64 



<*4S 

<*SS 



is a p, ^-staircase matrix. The characteristic of a staircase matrix is that if an 
element in the upper (lower) triangle is zero then all elements to the right (left) 
and above (below) are zero. Clearly, if A € S n , then p = 9; we say A is a 
p- staircase matrix. We are now ready for 



Theorem 6.6.3 If A € A/„ u NTN, it is a staircase matrix. 

Proof. The elements in the upper and lower triangles may be dealt with 
separately; we consider just the upper triangle. 

Suppose j > 1 and ay = 0. If k > j, then 

> 0, a,, = 0, a,, >0, aa> > 0 

imply a,* = 0; all elements to the right of a,; are afcx> zero. Now consider the 
terms {a,-,-, a, •_,•+!, . . . ,<*„) in the tth row. The first is positive; there is therefore 
a least index m, with i < m < n, such that a,, = 0 for j > m; call this index 
P.iP< > »• 

Now suppose j >i and a,, = 0. If k < » then 

> 0, an = 0, a„ > 0, a kj > 0 




fly <*1* 
a ii a >‘ 
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imply a k) = 0. Thus a,, = 0 and k < i implies a*j = 0. Thus j > p, implies 
a*,- = 0, i.e., j > p, implies j > p k ; p k < p,-. Thus p is a staircase sequence, and 
the upper triangle of A is a staircase. ■ 

Theorem 6.6.4 // A € M* is TN and 1 < p < n, then A(0(n);6(n)) £ 
AWp)Mp)) A(tf{py,tf{p)). Recall that 0{p) = {l,2,...,p} ( 0'(p) = {p + 

Proof. On account of Theorem 6.6.1 we may assume without loss of gener- 
ality that all the principal minors are positive, for if any were zero, then the in- 
equality would be satisfied trivially because then by Theorem 6.6.1, det(A) = 0. 
The theorem is true for n = 2 since 

A(0(2);0(2)) = a u a22 - a I2 u 21 < a,,^. 

We prow the theorem by induction, and assume that it holds for matrices of 
order n - 1 or less. We introduce the matrix B of Theorem 6.6.1: 

fey =A(0liii0uj), i,j = p+ l,...,n. 

= (A(9;0))-^ l A(0(ny,e(n)) 

= (D p )"~ p ~ 1 D n 

which we reverse to give 

D„ = /Wl/tV' - '- 

Since B(0 , |0 < | is of order n — p < n — 1, the inductive hypothesis applies to it: 

a I*';*') < b Pt i.y>iB(e\p + i);»'(p + i)) 

and thus 

Dn < by.^Bd / (p + 1);^(P + imw-*- 1 . (6.6.1) 

Applying Sylvester’s identity again, we have 

B{0\p + 1); 0 , (p + 1)) = 

where a = 0U (/(p + 1) = {1, 2, . . . ,p, p + 2, . . . , n} which when combined with 
(6.6.1) and fe ptl .,,»i = D ptl gives 

Dn < 

< Dy,\A{a\a)/D v ' (6 • 6 • 2, 

Now we use the inductive hypothesis again to give 

A{a:«) < D p A{#(p+ l);^(i>+ 1)) 
which, when combined with (6.6.2) gives 

A(*(»);*(n)) < 1 );%+ l);^(p+ 1)) 



which shows that the result holds for matrices of order n. ■ 
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Corollary 6.6.2 //A € M„ is TN then 

Dp < 011022 . - • Oyp- 1 < P < «• 

Theorem 6.6.4 is expressed as a result concerning principal minors of a TN 
matrix A, but since any square matrix taken from a subset of rows and columns 
of such an A is also TN we can state 

Corollary 6.6.3 // A € Af n is TN, 0,y € Q,. n , and 0,-y € 9(p) = t) (i.e., 

0 r 1 q Zp), 

A(0ue , ^ue , )<A{0^)A^ l e , j. 

Similarly, if € Q q .», and 0,i G t>'(p ) = f? (i.e., 0,,?, > p + 1), then 

A{OU0-,OUi)<A{0-,y)A(O:O). 

See Ando (1987) (4| for generalisations of this result. 

Theorem 6.6.5 Suppose A G A/ m n is TN. If A has p linearly dependent rows, 
labelled by a € Q p , n teith a I = l,a r = m, of uhich the first p — 1, labelled o\o v , 
and the last p— 1, labelled a\ai, are linearly independent fit.), then A has rank 
P-1- 

Proof. Clearly, the rank of A is at least p — 1; we show that it is not greater 
than p — 1, i.e., it is exactly p — 1. 

The linearly dependent rows are specified by a = (aj, q 2 , . . . , a,,}. If p > n, 
then rank (A) < n < p, so that rank (A) = p - 1. Take p < n. The row Op 
may be expressed in terms of rows ai,a 2 i - - • > Q p-i : 

p-i 

•W.r = 53 c * fl a, j 3 = 1. 2,-. - ,n, (6.6.3) 

1-1 

and since rows a\aj are C\ ^ 0. Since rws o\Oj| are l.L, there is tfl € 
Q p -\ %n such that > 0. On substituting for o 0#0 from (6.6.3) we 

find 

A(a\« l ;0°> = (-rc I A(a\«,;0 o ). 

Therefore A(a\ai;fP) ^ 0, but this minor is non-negative and therefore it is 
strictly positive; therefore {-) p c x > 0. 

Now suppose q € a' then < q < a r +i for some index r satisfying 1 $r< 
p — 1. If p € Q Pt n then, on substituting for a 0 , ti , as before, we find 

^i)U^) = (-^'ciAU^aJUq'd). (6.6.4) 

The inequality (-) p c\ > 0 implies that the minors on either side of (6.6.4) have 
opposite signs. But both are non- negative so that both are zero. Since p is an 
arbitrary member of Qp.n, this means that any row q Q d may be expressed as 
a linear combination of the rows a\oi, or equivalently of a\a p . Thus the rank 
of A is p — 1. ■ 

We now prove a corollary' of this result, but since its truth is not immediately 
clear, we state it as 




